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Matrix partition of U2C +D r 
with position partials. 
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9.1-3 


10.2-4 


Acceleration of satellite due to drag 3.1-2 
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Surface density acceleration 3.3-28 
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Weighted RMS of previous outer iteration 
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10.3-2 
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11.4-1 
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I Identity matrix 9.1-3 
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nutation calculation 
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11.4-1 
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8.3-29 
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8.3-27 
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Direction cosine (measurement type) 

6.1-7 

®d 

Mass of .the disturbing body for third 
body perturbations 

8.4-1 

“i 

Computed equivalent of the i th measure* 
ment (see C ^ and 

10.2-2 

“s * 

Mass of the satellite 

8.5-4 
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Number of observations in £ 

10.1-1 

N' 

Maximum degree coefficient unaffected 
by. the surface density layer 

8.3-27 
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North baseline unit vector in the S.2-1 

topocentric horizon coordinate system 

i 

Nj.- Residuals contributing to the bias 11.3-1 

computation 

n Direction cosine (measurement type) * 6.1-7 

n Number of residuals 11.3-1 

n^ Number of electronic biases 11.3-1 

n Surface index of refraction 7.5-5 

s 

0^ The i** 1 observed measurement value 2.2-1 

F Vector of parameters to be determined 2.2-1 

( ) Legendre polynomial 6.3-2 

P Solar radiation pressure in the vicinity 8.5-4 

of the Earth 

p (x) Joint probability density function x 10.1-1 

* 

p (z) Joint probability density function for £ 10.1-1 

p (xjz) Joint conditional probability density 10. 1-1 

function for x, given that z has 
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Symbols 
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p (z|x) Joint conditional probability density 
function for z_ give.* that x has 
occurred 


q Parallactic angle in radians 

R e Mean earth radius 

Third body disturbing potential 

R^ Distance from center of mass of the 

earth to the center of mass of the 
disturbing body 


Rg(t) Range vector from the center of the 

earth to the ground station r.t time t 


R i Residual value (dm^) 

Unit vector from center of mass of the 
earth to the center of mass of the 
disturbing body 

R,(t) Range sum measurement at time t 

ft 

R^(t) Range vector from the center of the 

earth to the relay satellite at time t 


R 2 (t) 


Range vector from the center of the 
earth to the tracked satellite at time 
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Symbols Description First Used 

Down-link range from the relay satellite 6.4-3 
to the ground 


R lu Dp-link range from the ground to the relay 6.4-1 

satellite 




R 2d 

Relay satellite - tracked satellite range 

6.4-1 

R 2u 

Tracked satellite - relay satellite 
range 

6.4-3 

R s 

Time derivative of R $ 

6.4-8 

»lu 

Time derivative of R^ u 

6.4-8 

R 2d 

Time derivative of 

6.4-8 

r 

Distance from the point of interest to 
the center of mass of the earth 

8.3-18 

r 

Distance from center of mass of the 
earth to satellite 

8.8-1 

r 

Geocentric satellite position vector 

5.1-10 

r 

True of date position vector of the 
satellite 

8.7-29 


?£ True of date position vector of third 

body for third body gravitational effects 
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Page 

Description First Used 

Geocentric position vector of a tracking 2.2-4 
station 


a 

r 


A 



s 


s 


Unit vector from center of mass of the 8.8-1 
Earth to the satellite 

True of date unit vector pointing to ' 8.5-4 

the Sun 

- . . . * 

The cosine of the enclosed angle between 8.4-1 
r and r^ 

Surface of the Earth .8.3-18 


S 1 

The first sum carry along by the 
integrator 

9.3-1 

S 2 

The second sums carry along by the 

9.3-1 


integrator , . .. 

h . _ -t. 

S nm 

Gravitational harmonic coefficient 
of degree n , order m 

8.2-2 

**nm 

The sine coefficient of surface density 
constraint equations . ... w 

8.3-25 

, 2 

• * * ; " 

Sample variance 

11.3-1 

T 

A sample layer distributed on the 
surface of the Earth 

8.3-25 
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T Exospheric temperature 

T Exospheric temperature 


Page 

First Used 

8.7- 15 

8.7- J 
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1 > Sl 
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* 


u 

u 
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2C 


A 

V 


Nighttime minimum global exospheric 
temperature for a given day 

Average nighttime minimum global 
exospheric temperature for a given 
period 

Geopotential field of the Earth 

Spherical harmonics part of total 
earth potential 

Matrix containing the second partial 
derivatives of the gravitational 
potentials with respect to the true 
of date position coordinates 

Central angle between the satellite 
vector and a vector pointing toward 
the ascending node of the orbit 

Unit vector in the direction of p" 

A 

Covariance matrix of x 

Unit local vertical at the station 
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Page 
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V A 

a priori covariance matrix associated 

10.2-1 


with x^; same as 


V a 

Matrix partition of V^; a priori co- 

i 

10.2-3 


variance matrix associated with a 


V k 

Matrix partition of V^; a priori co- 

10.2-3 


variance matrix associated with k 


V r 

Matrix partition of V ft associated 
vith the r** 1 arc 

10.2-6 

W 

Weighting matrix for observations; 
-1 

same as E, 
z 

10.2-1 

W 

Total potential of the Earth 

8.3-.18 

X 

C 'ordina't : system direction : 

2.1-3 


a) Direction in the equatorial plane 
pointing toward the Greenwich 
meridian (Earth-fixed system) 


b) In the direction of the true equinox 
of date at o^o of the epoch day 
(inertial system) 

c) In the direction of the true equinox 
of date (true of date system) 
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X(t+At) Position partial at time t 


9.3-1 


X(t+At) Velocity partial at time t 

X The X angle of the satellite 

(measurement type) 

* 

X g Earth-fixed position component 

X. True of date position component 

X m Matrix containing the variational 

partials 

Xj^ Inertial cartesian position coordinates 

of the relay satellite 

X 2i Inertial. cartesian position coordinates 

of the tracked satellite 

X u Time derivative of X^ 

Time derivative of Xj^ 

x True of date X position component 

of the satellite 

x Rotation angle for polar motion 


9.3-2 


6.1-7 


3.4-1 


3.4-1 


8 . 2-6 


6.4-6 


6.4-6 


6.4-8 


6.4-8 


2.2-4 
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(n) 
*A 


Description 

Vector of M parameters 

The ’’best” estimate of x 

t* h A 

The n approximation to x 

The a priori estimate of x 

The vector describing the true of 
date position and velocity of the 
satellite 


Y Coordinate system direction 

(associated with the X and Z 
directions) 


Y Partition of X m ; a matrix containing 
3r 

• ^ 

Y Partition of X m ; a matrix containing 
3 ? 

•• *• 

Y Matrix containing 3jF; same as 

matrix F 

* 

Y a The Y angle of the satellite 

(measurement type) 

Y e Earth-fixed position component 
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First Used 
10 . 1-1 

10 . 1-2 

10 . 1-2 

10 . 1-2 

2 . 2-4 

2 . 1 - 3 

9 . 1 - 3 

9 . 1 - 3 

9 . 1 - 3 

6 . 1 - 7 
3 . 4-1 
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True of date position component 

True of date Y position 
component of the satellite 

Rotation angle for polar motion 

* 

Direction of the spin axis of the Earth 
for Z direction of coordinate systems. 
(Taken at o^o of epoch day for inertial 
coordinate system.) Compare X 

The zenith baseline unit vector in the 
topocentric horizon coordinate system 

Earth- fixed component; same as z 

Observed zenith angle 

True of date Z position coordinate 
of the satellite 

A precession angle 

A vector of N independent observations 


Page 

First Used 

3.4- 1 

2.2- 4 

5.4- 5 
2 . 1-2 

5.2- 1 

5.1- 5 

7.5- 1 

2.2- 4 

3.1- 1 

10 . 1 - 1 
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o Topocentric right ascension of the 

satellite (measurement type) 

✓ 

a" Observed declination of the satellite 

o The set of parameter not affecting 

the dynamics of satellite motion 

£ The set. of parameters affecting the 

dynamics of satellite motion 

y Parameter of differential corrections 

for epoch element and force model parame 
ter errors 




Cowell integration scheme coefficients 


AE^ Area of the surface density block 

A* Correction to measurement of direction 

cosine A 


Am Correction to measurement of direction 

cosine m 


AR Differential refraction 
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6 . 1 - 5 
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3 . 3 - 18 

7 . 5 - 5 
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Page 
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AT 

09 

Geomagnet ic beating correction to T w 

8.7-9 

At 

Measurement timing bias 

6.0-2 

^ i 

4t i« 

Transit time for the range R ly 

6.4-3 

4t ia 

Transit time for the range R^. 

6.4-3 

6t 2d 

Transit time for the range R 2 ^ 

6.4-3 

At 2u 

Transit time for the range R 2 U 

6.4-3 

ix a 

Correction to measured X angle 

7.6-1 

AY « 

Correction to measured Y angle 

7.6-1 

Act 

Equation of the equinoxes 

3.5-2 

Aa 

Right ascension measurement correction 

7.4-1 

A6 

Declination measurement correction 

7.4-1 

Ae 
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SECTION 1.0 
THE GEODYN PROGRAM 


GEODYN was written for GSFC by WOLF in 1971 and has 
been operational since January of 1972. A merger of the 
Multi -Arc NONAME program and the GE03TAR program, GEODYN is 
greatly improved in overall capability, accuracy, and versa- 
tility over its parent programs. 

GEODYN is one of the most widely used orbit and geo- 
detic parameter estimation programs in the world. It is 
currently operational at GSFC on the IBM 360 *95, '91, and 
*75; at Onio State University on the IBM 370/155; and will 
shortly be operational at Wallops Island on the GE 635. 
Additionally the GEODYN parent program. Multi -Arc NONAME 
is operational at the Goddard Institute for Space Studies 
in New York on an IBM 360/95 and at the Institut fur Physik 
and Plasmaphysik, Garching, West Germany on an IBM 360/51. 

GEODYN has been used for 

• determination of definitive orbits 

• tracking instrument calibration 

• satellite operational predictions 

• geodetic parameter estimation 

and many other items relating to applied research in 
satellite geodesy using virtually all types of satellite 
tracking data. 


1.0-1 


SECTION 2.0 

THE ORBIT AND GEODETIC PARAMETER ESTIMATION PROBLEM 


The purpose of this section is to provide an under- 
standing of the relationship between the various elements 
in the solution of the orbit and geodetic parameter esti- 
mation problem. As such, it is a general statement of 
the problem and serves to coordinate the detailed solutions 
to each element in the problem presented in the sections 
which follow. 

The problem is divided into two parts: 

• the orbit prediction problem, and 

• the parameter estimation problem. 

The solution to the first of these problems corresponds to 
GEODYN's orbit generation mode. The solution to the 
latter corresponds to GEODYN's data reduction mode and 
of course is based on the solution to the former. 

• 

The reader should note that there are two key choices 
which dramatically affect the GEODYN solution structure: 

• Cowell's method for integrating the orbit, and 

• a Bayesian least squares statistical estimation 
procedure for the parameter estimation problem. 


2.1 


The Orbit Prediction Problem 


There are a number of approaches to orbit pre- 
diction. The GEODYN approach is to use Cowell’s method, 
which is the direct numerical integration of the satellite 
equations of motion in rectangular coordinates. The 
initial conditions for these differential equations are 
the epoch position and velocity; the accelerations of 
the satellite must be evalu^Led. 

The acceleration 7 oducing forces which are cur- 
rently modelled m GEODYN are the effects of 

o the geopotential, 

o surface densities, 

o the luni-solar potentials, 

o planetary potentials of Venus, Mars, Jupiter 
and Saturn, 

o Radiation pressure, 

♦ 

o earth tidal potential, and 

0 atmospheric drag 

Perhaps the most outstanding common feature of these forces 
is that they are functions of the position of the satellite 
relative to the Earth, Spn, Moon, or. Planets and of the Sun 
and Moon relative to the Earth. Only atmospheric drag is a 
function of any additional quantity,* specifically, the rela- 
tive velocity of the satellite with respect to the atmosphere 


•Not to be confused with the ’’fixed" parameters in the models 
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The accurate evaluation of. the acceleration of 
a satellite therefore involves the solution to two 

v\ ^ 4 4* mvi 4> rt v> U 1 a tw • 
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o the accurate modeling of each force on the 
satellite - Earth - Sun - Moon - Planet 
relationship, and 

o The precise modeling of the motions of the 
Earth, Sun, Moon, and Planets. 

The specific details for each model in these solutions 
are given elsewhere in Sections 3, 4, and 8. The 
question of how these models fit together is in effect 
the question of appropriate coordinate systems. 

The key factor in the selection of coordinate 
systems for the satellite orbit prediction problem is 
the motion of the Earth. For the purposes of GEODYN, 
this motion consists of: 

o precession and nutation, and 

o rotation. 

We are considering here the motion of the solid body o. 
the Earth, as versus the slippage in the Earth's crust 
(polar motion) which just affects the position of the 
observer. 

The precession and nutation define the variation 
in 


o 


the direction of the spin axis of the Earth 
(♦ Z), and 



o the direction of the true equinox of date 
(+ X). 

These directions define Ihe (geocentric) true of date 
coordinate system. 

The rotation rate of the Earth is the time rate 
of change of the Greenwich hour angle 9^ between the 
Greenwich meridian and the true equinox of date. Thus 
the Earth-fixed system differs from the true of date 
system. according to the rotation angle 8 . 

o 

The equations of motion for the satellite must be 
integrated in an inertial coordinate system. The GEODYN 
inertial system is defined as the true of date system 
corresponding to o !?0 of a reference epoch. 

The coordinate systems in which the accelerations 
due to each physical effect are evaluated should be 
noted. The geopotential effects are evaluated in the 
Earth- fixed system, and then transformed to true of 
date to be combined with the other effects. The others 
are evaluated in the true of date system. The total 
acceleration is then transformed to the reference inertial 
system for use in the integration procedure. 

The integration procedure used in GEODYN is a 
predictor- corrector type with a fixed time step. There 
is an optional variable step procedure. As the integration 
algorithms used provide for output on an even step, an 
interpolation procedure is required. 


2.2 


The Parameter Estimation Problem 


Let us consider the relationships between the 
observations CL, their corresponding computed values 
and P, the vector of parameters to be determined. These 
relationships are given by- 



CD 


where 


i 


denotes the i** 1 observation or association 
with it. 


dPj is the correction to the j** 1 parameter, and 

dO. is the error of observation associated with 
1 

the i * observation. 


The basic problem of parameter estimation is to determine 
a solution to these equations. 

The role of data preprocessing is quite apparent 
from these equations. First, the observation and its 
computed equivalent must be in a common time and spatial 
reference system. Second, there are certain physical 
effects such as atmospheric refraction which do not 
significantly vary by any likely change in the parameters 
represented by F. 


These computations and corrections may equally 
well be applied to the observations as to their computed 


values. Furthermore, the relationship between the computed 
value and the model parameters F is, in general, nonlinear, 
and hence the computed values may have to be evaluated 
several times in the estimation procedure. Thus a con- 
siderable increase in computational efficiency may be 
attained by applying these computations and corrections 
to the observations; i.e., to preprocess the data. 

The preprocessed observations used by GEODYN are 
directly related to the position ^nd/or velocity of the 
satellite relative to the observer at the given observa- 
tion time. These relationships arc geometric, hence 
computed equivalents for these observations are obtained 
by ?nplying these geometric relationships to the computed 
values for the positions and velocities of the satellite 
and the observer at the desired time. 

Associated with each measurement from each ob- 
serving station is a (known) statistical uncertainty. 

This uncertainty is a statistical property of the noise 
on the observations. This uncertainty is the reason 
a statistical estimation procedure is required for the 
GEODYN parameter determination. 

It should be noted that dO^, the measurement 
error, is not the same as the noise on the observations. 

The dO^ account for all of the discrepancy which 

is not accounted for by the corrections to the parameters 
c*F. These dO^ represent both 

t the contribution from the noise on the 
observation, and 

• the incompleteness of the mathematical model 
represented by the parameters F; 


By this last we mean either that the parameter set 
being determined is insufficient or that the functional 
form of the model is inadequate. 

GEODYN has two different ways of dealing with 
these errors of observation: 

1. The measurement model includes both a 
constant bias and a timing bias which may 
be determined. 

2. There is an automatic editing procedure 
to delete bad (statistically unlikely) 
measurements . 

The nature of the parameters to be determined has 
a significant effect on the functional structure of the 
solution. In GEODYN, these parameters are: 

• the position and velocity of the satellite 
at epoch. These are the initial conditions 
for the equations of motion. 

• force model parameters. These affect the 
motion of the satellite. 

• station positions and biases for station 
measurement types. These do not affect the 
motion of the satellite'. 

Thus, the parameters to be determined are implicitly 
partitioned into a set a, which are not concerned with 
the dynamics of the satellite motion and a set F which 
are. 
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The computed value for each observation 0^ is a 


function of 


"ob *^e E art h- fixed position vector o r the 
station, and 


x t true of date position and velocity vector 


of the satellite (x,y ,z ,x,y ,z} 


at the desired observation time. When measurement biases 
are used, is also a function of F, the biases associated 
with the particular station measurement type. 


; I 


1 ; 


Let us consider the effect of the given partitioning 
on the required partial derivatives in the observational 


■ ( 

equations. 
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The partial derivatives are called the variational 
partials. While the other partial derivatives on the 
Tight-hand side of the equations above are computed from 
the measurement model at the given time, the variational 
ptitials must be obtained by integrating the variational 
equations. As will be shown in Section 8, these equations 
are similar to the equations of motion. 


ii t 
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1 


i 


Tae need for the above mentioned variational 
partials obviously has a dramatic effect on any solution 
to the observational equations. In addition tc integrat- 
ing the equations of motion to -enerate an orbit, the 
solution requires that the variational equations be 
integrated. 

We have heretofore discussed the elements of the 
observational equations; we shall now discuss the solution 
of these equations; i.e., the statistical estimation 
scheme. 


There are a number of estimation schemes that 
can be used. The method used in GEODYN is a batch 
scheme that uses all observations simultaneously to 
estimate the parameter set. The alternative would be 
a sequential scheme that uses the observations se- 
quentially to calculate an updated set of parameters 
from each additional observation. Although batch and 
sequential schemes are essentially equivalent, practical 
numerical problems often occur with sequential schemes, 
especially when processing highly accurate observations. 
Therefore, a batch scheme was chosen. 

The particular method selected for GEODYN is a 
partitioned Bayesian least squares method as detailed 
in Section 10. A Bayesian method is selected because 
such a scheme utilizes meaningful a priori information. 
The partitioning is such that the arrays which must be 
simultaneously in core are arrays associated with 
parameters common to all satellite arcs, and ^rvays 
pertaining to the arc being processed. Its purpose is 
to dramatically reduce the core storage requirements ot 
the program without any significant cost in computation 
time . 
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There is an interesting aside related to the use 
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information for the parameters guarantees that the esti- 


mation procedure will mechanically operate (but not 
necessarily converge). The user must ensure that nis 
data contains information relating to the parameters 
he wishes determined. 
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SECTION 3.0 

THE MOTION OF THE EARTH AND RtLAiED COORDINATE SYSTEMS 


ftSPROKLCteaj,. 
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The major faccor in satellite dynamics is the 
gravitational attraction of the Earth. Because of the 
(usual) closeness of the satellite and its primary, 
the Earth cannot be considered a point mass, and hence 
any model for the dynamics must contain at least an 
implicit mass distribution. The concern of this section 
is the motion of this mass distribution and its relation 
to coordinate systems. 

We will first consider the meaning of this motion 
of the Earth in terms of the requisite coordinate systems 
for the orbit prediction pioblem. 

The choice of appropriate coordinate systems is 
controlled by several factors: 

• In the case of a satellite moving in the 
Earth's gravitational field, the most 
suitable reference system for orbit com- 
putation is a system with its origin at 
the Earth's center of mass, referred to 
as a geocentric reference system. 

• The satellite equations of motion must be 
integrated in an inertial coordinate system. 

, * 

• The Earth is rotating at a rate 0^, which is 

the time rate of change of the Greenwich hour 
angle. This angle is the hour angle of the 
true equinox of date with respect tc the 
Greenwich meridian as measured in the equatorial 
plane . 
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» The. Earth both precesses and nutates, thus 
changing the oiretliuns of both the Earth s 
spin axis ~nd the true equinox of date in 
inertial space. 

The motions of the Earth referred to here are of course 
those of the "solid body" of the Earth, the motion of 
the primary mass distribution. The slippage of the 
Earth's crust is considered elsewhere in Section 5.2 
(polar motion). 


3.1 The True of Date Coordinate System 

Let us consider that at any given time, the spin 
axis of the Earth (* Z) and the direction of the true 
equinox of date (+ X) may be used to define a right-handed 
geocentric coordinate system. This system is known as 
the true of date coordinate system. The coordinate 
systems of GEODYN will be defined in terms of this system. 

REPRODUCIBILITY OF THE 

3 . 2 The Inertial Coordinate System URiGINAL PAGE IS POOR 

The inertial coordinate system of GEODYN is the 
true of date coordinate system defined at ($0 of the 
reference day for each satellite. This is the system in 

which Che satellite equations of motion are integrated. 

* 

This is a right-handed, Cartesian, geocentric 
coordinate system with the X axis directed along the true 
equinox of 0.0 of the reference day and with the Z axis 
directed along the Earth's spin axis toward north at the same 
time. The Y axis is of ~ourse defined so that the co- 
ordinate system Is orthogonal. 
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It should be noted that the inertial system differs 
from the true of date system by the variation in time of 
the directions of the Earth's spin axis and the true 
equinox of date. This variation is described by the 
effects of precession and nutation. 


3.3 The Earth- fixed Coordinate System 

The Earth-fixed coordinate system is geocentric, 
with the Z axis pointing north along the axis of rotation 
and with the X axis in the equatorial plane pointing 
toward the Greenwich meridian. The system is orthogonal 
and right-handed; thus the Y axis is automatically defined. 

This system is rotating with respect to the true 
of date coordinate system. The Z axis, the spin axis of 
the Earth, is common to both systems. The rotation rate 
is equal to the Earth's angular velocity. Consequently, 
the hour angle 9^ of the true equinox of date with respect 
to the Greenwich meridian (measured westward in the equa- 

o 

torial plane) is changing at a rate 6 equal to the angular 

8 

velocity of the Earth. 



3.4 Transformation Between Earth-fixed and True of 
bate Coordinates ' 

, ♦ 

The transformation between Earth- fixed and true of 
date coordinates is a simple rotation. The Z axis is 
common to both systems. The angle between X^, the true 
of date X component vector, and X A , the Earth-fixed 
component vector, is 9^, the Greenwich hour angle. The 
Y component vectors are similarly related. These trans- 
formations for X 0 , Y e , X^, Y^ which ere accomplished in 
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GEODYN by the functions XEFIX, YEFIX, XINERT, and YINERT GRHRAN 
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The transformation of velocities requires taking 
into account the rotational velocity, 0 g , of the Earth- 
fixed system with respect to the true of date reference 
frame. The following relationships should be noted: 
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The velocity transformations are then 
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The brackets denote the part of each transform which is 
a transformation identical to its coordinate equivalent. 


These same transformations are used in the 
transformation of partial derivatives from the Earth- 
fixed system to true of date. For the k measurement, 
C^, the partial derivative transformations are ex- 
plicitly: 
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The brackets have the same meaning as before. 

These above transforms are used or computed 
using the functions XEFIX, YEFIX, XINERT, or YINERT 
in three GEODYN subroutines: GRHRAN, OBSDQT, and 

PREDCT . 


3.5 Computation of 6 
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The computation of the Greenwich hour angle is quite GRHRAN 
important because it provides the orientation of the Earth F 
relative to the true of date system. The additional affects; 
i.e., to transform from true of date to inertial, of pre- 
cession and nutation are sufficiently small that early orbit 
analysis programs neglected them. Thus, this angle is the.,, 
major variable in relating the Earth- fixed system to the 
inertial reference frame in which the satellite equations 
of motion are integrated. -v-: - 
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The evaluation, of 9^ is disr.is«.r-d in detail in GRHRAM 

the Explanatory Supplement , Reference 1. 6^ is computed F 

in subroutines GRHRAN and F from the expression: 


B * e + At. 0. + At- 0- + Aa 
g 11 2 2 


Atj is the integer number of days since 
January 0.0 UT of the reference year, 

A ^2 is the fractional UT part of a day for 
the time of interest, 

6^ is the Greenwich hour angle on 

0 January 0.0 UT of the reference year, 


8^ is the mean advance of the Greenwich 
hour angle per mean solar day. 


0£ is the mean daily rate of advance of 
Greenwich hour angle and 

Aa is the equation of equinoxes (nutation in 
right ascension) . 


The initial 6 is obtained from a table of 

g 

values containing the° Greenwich hour angle 
on January 0.0 for each year. This table is in 
Common Block CGEOS and is accessed in JANTHG. 


JANTHG 
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The equation of equinoxes, Aa, is obtained from 
subroutine EPHEM, which calculates the quantity from 
the ephemeris tape data according to the Everett fifth- 
order interpolation scheme. 
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Precession and Nutation 
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The inertial coordinate system of GEODYN, in NUTATE 

vhich the equations of motion arc integrated, is de- PRECES 

fined by the true equator and equinox of date for REFCOR 

0^0 of the reference day. However, the Earth-fixed 
coordinate system is related to the true equator 
and equinox of date at any given instant. Thus, it is 
necessary to consider the effects which change the 
orientation in space of the equatorial plane and the 
ecliptic plane. 

These phenomena are 

9 the combined gravitational effect of the 
moon and the sun on the Earth's equatorial 
bulge, and 

• the effect of the gravitational pulls of 
the various planets on the Earth's orbit. 

The first of these affects the orientation of the 
equatorial plane; the second affects the orientation 
of the ecliptic plane. Both affect the relationship 
between the inertial and Earth-fixed reference systems 
of GEODYN. 

The effect of these phenomena is to cause pre- 
cession and nutation, both for the spin axis of the 
Earth and for the ecliptic pole. This precession and 
nutation provides the relationship between the inertial 
system defined by the true equator and equinox of the 
reference date and the "instantaneous" inertial system de- 
fined by the true equator and equinox of date at any 
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given instant. Let us consider the effect of each of 
these phenomena in greater detail. 

The luni-solar effects cause the Earth's axis 
of rotation to precess and nutate about the ecliptic 
pole. This precession will not affect the angle be- 
tween the equatorial plane and the ecliptic (the 
"obliquity of the ecliptic""' but will affect the 
position of the equinox in the ecliptic plane. Thus 
the effect of luni-solar precession is entirely in 
celestial longitude. The nutation will affect both, 
consequently we have nutation in longitude and nuta- 
tion in cb’ iqui ty. 

The effect of the planets on the Earth’s orbit 
will cause both secular and periodic deviations. 
However, the ecliptic is defined to be the mean plane 
of the Earth’s orbit. Periodic effects are not con- 
sidered to be a change in the orientation of the 
ecliptic; they are considered to be a perturbation 
of the Earth's celestial latitude. (See Reference 1.) 

The secular effect of the planets on the 
ecliptic plane is separated into two parts: planetary 

precession and a secular change in obliquity. The 
effect of planetary precession is entirely in right 
ascension. 

In summary, the secular effects on the orienta- 
tions of the equatorial plane are: 
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PRECES 

REFCOR 



luni-sclar precession, 


• planetary precession, and 

• a secular change in obliquity. 

As is the convention, all of these secular effects are 
considered under the heading, ’’precession." The 
periodic effects are 

• nutation in longitude , and 

• nutation in obliquity. 

In terms of the GEODYN system, subroutine PRECES 
determines the secular effects; i.e., the rotation 
matrix which will transform coordinates from the mean 
equator and equinox of date to the mean equator and 
equinox of 1950.0. 

Subroutine NUTATE determines the rotation matrix 
to transform from true equator and equinox of date to 
mean equator and equinox of date. This accounts for 
the periodic effects. 

GEODYN has two different routines for transform- 
ing from one epoch to another. These are EQUATR and 
REFCOR. EQUATR will take either mean or true coordinate 
input and will output in either mean or true coordinates, 
REFCOR will take only true coordinate input and will output 
only true coordinates. The same general algorithm is used 
in both: 


EQN 

EQUATR 

NUTATE 

PRECES 

REFfOR 


PRECcS 


NUTATE 


EQUATR 

REFCOR 


- 4 * 


o Rotate from true to mean equator and 
equinox of input date if required. 

o Rotate from mean of input date to mean 
of 1950.0. 

o Rotate from mean of 1950.0 to mean 
of output date. 

o Rotate from mean to true of output date 
if required. 

All of these rotations are of course done with rotation 
matrices. 

Subroutine REFCOR will transform between any 
time of day and 0*?0 on a given reference day. It performs 
this transform by interpolating linearly between the ro- 
tation matrices for the day of the input and that day plus 
one . 


3.6.1 Precession 

The precession of coordinates from the mean 
equator and equinox of one epoch tQ to the mean equator 
and equinox of is accomplished very simpiy. Ex- 
amine Figure 1 and consider a position described by the 
vector X in the X^ ,X 2 ,Xj coordinate system which is 


i-’ r 
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defined by tbe mean equator and equinox of tg. Like- 
wise, consider the same position as described by the 
vector 7 in the system defined by the mean 

equator and equinox of t^ . The expression relating 
these vectors, 


7 = R 3 (- 2 ) R 2 (0) R 3 (-O 5C, (1) 

follows directly from inspection of Figure 1. 

It snvUld be observed that 90° - s is the 
right ascension of the ascending node of the equator 
of epoch tg reckoned from the equinox of tg, 90° - z 

is the right ascension of the node reckoned from tne 

equinox of t^ and 9 is tne inclination of the equator 
of to the epoch of tg. 

Numerical expressions for these rotation angles 
2,6, C were derived by Simon Newcomb, based partly upon 
theoretical considerations but primarily upon actual 
observation. (See References for the derivations.) 

The formulae used in GEODYN are relative to an initial 
epoch of 1950.0; 

C » ?30S 953 204 65 x 10' 6 d ♦ ?109 749 2 x 10" 14 d 2 (2) 

♦ ?178 097 x I0‘ 20 d 3 

2 - ?305 953 204 65 x 10' 6 d + *397 204 9 x 10' 14 d 2 (3) 

+ *191 031 x 10* 20 d 3 
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n )Lte: 

4 . 


1 

i 

l 




I 
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3.6.2 Nutation 


NUTATE 


The nutation of coordinates between mean and 
true equator and equinox of date is readily accomplished 
using rotation matrices. Examine Figure 1 and consider 
a position described by the vector X in the X- .X ,X 

A • J 

system which is described by the mean equator and equi- 
nox of date. Likewise, consider the sane position as 
described by the vector Z in the s y stem de- 

fined by the true equator and equinox o£ date. The 
expression relating these vectors, 


f 



Z . Rj C-e T ) Rj (-4*) Rj (cj X, (1) 


follows directly from inspection of Figure 1. 

The definition cf these angles are: 
e.^ - true obliquity of date 
e *- mean obliquity of date 
Af - nutation in longitude 

Note that - e ffl is the nutation in obliquity. 

* 

The remaining problem is to compute the nutations 
iri longitude- and obliquity. The algorithm used in 
GEODYN was developed by Woolard and is coded in sub- 
routine EQN. 


NUTATE 

EQN 


5.6-3 





NUTATION 


7 . x, 

■ -5 






*M \ 

Tr J# E^oWOt^S 


8 Mean Ooliquity of Date 

« T « True Obliquity of Date 

Y h 8 Direction of Mean Equinox of Date 

Y t ■ Direction of Time Equinox of Date 


Figure It Rotation Between Mean Equator & Equinox of Date 

and 

True Equator & Equinox of Date 


Woolard's solution as it appears m references i 
through 4 is reproduced in Tables la, lb, and lc. The 
periodic terms have been rearranged in descending order 
of magnitude. The subprogram EQN computes the nutation 
in longitude and obliquity by using th^ algorithm in 
Tables 2a, 2b, and 2c. Ir Table 2a the angular units 
of the fundamental arguments have been changed to 
radians and the time units have been changed to days. 
Tables 2b and 2c are identical to Tables lb and lc 
..ften neglecting ail periodic terms with coefficients 
less than 7001 and all secular portions of the co- 
efficient which are less than 7001. The expressions 
for true obliquity of date and nutation in right as- 
cension appear in Table 2d. 

The definitions of the variables used in these 
solutions and additional notation are as follows: 

J * Julian Ephemeris Date of desired calculation 

J ■ 241 S020.5 (Julian Ephemeris Date corresponding 

to 1900 January 0.5 Epnemeris Time) 

T » (J-J 0 )/36525 ■ Julian ephemeris centuries of 

3652S Ephemeris Days elapsed from J Q to J 

d * J~J 0 “ Ephemeris Days elapsed from J to J Q 


tyjy 
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COORDINATE SYSTEM: Geocentric, eel ip tic ami 

mean equinox of date: 

g « mean anomaly Moon 

g' * mean anomaly - Sun 

c *» mean angular distance of the Moon from its 
ascending node 

D * mean elongation of the Moon from the Sun 

ft * longitude of the mean ascending node of the 

Moon's orbit 

e u « mean obliquity of date 

C 

e.j. ■ true obliquity of date 

Ae « nutation in obliquity 

Ai| > ■ nutation in longitude 

Aa ■ nutation in right ascension 
(equation of the equinoxes) 


I'All 

L, yi« 


f 




REPRODUCIBILITY of THb 
ORIGINAL PACE IS POOR 


TABLE la FUNDAMENTAL ARGUMENTS 


g * 296 

°06"16VS9 

+ 

1325 r 198° 

50' 56 

V 7 9 T 

4 

33V09 

T 2 + 

VO 

518 

I 3 

g ' *3358 

° 28 ' 33V 00 

+ 


99 r 359° 

02*59 

V10 T 

- 

V59 

T 2 - 

"0120 

T -5 

F - 11 

0 15 ‘ 03V 20 

+ 

1342 r 82° 

01*30 

VS4 T 

- 

11V56 

T 2 - 

V0012 

T 3 

D = 350 

°44 ' 14V95 

+ 

12 

36 r 307° 

06*51 

V18 T 

- 

5 V 1 7 

T 2 - 

V0068 

T 3 

ft » 259 

°1C ' 59V 79 

- 


5 r 134° 

08 ' 31 

V23 T 


7V48 

T 2 + 

*.’0080 

T 3 

V 23 

°27‘ 08V26 

- 



46 

V845T 

- 

V00S9T 2 + 

'.*0080 

T 3 



TABLE lb 

NUTATION IN OBLIQUITY 















Series 

No. 

Ae •• + 

(+0V00091 

T 

4 

9V2100) 

COS 

( 





♦ 

ft) 

1 

4 

(-0V00029 

T 

4 

0.5522) 

cos 

( 


4 

2F 

- 2D 

4 

2ft) 

2 

4 

(+0.00004 

T 

- 

0.0904) 

cos 

( 





4 

2ft) 

3 

4 

(-0.00005 

T 

4 

0.0884) 

cos 

( 


4 

2F 


4 

2Q) 

4 

4 

(-0.00006 

T 

4 

0.0216) 

cos 

( 

4 

g* ♦ 

2F 

- 2D 

4 

2ft) 

5 




4 

0.0183 

cos 

( 


4 

2F 


4 

O) 

6 

+ 

(-0.00001 

T 

♦ 

0.0113) 

CCS 

(+ g 


♦ 

2F 


4 

2ft) 

7 

4 

(+0.00003 

T 

- 

0.0093) 

cos 

( 

- 

g’ + 

2F 

- 2D 

4 

• ') 

8 




- 

0.0066 

cos 

( 


4 

2F 

- 2D 

4 

ft) 

9 


. 


- 

0.0050 

cos 

(* g 


4 

2F 


4 

2ft) 

10 




- 

0.0031 

cos 

( + g 





4 

fti 

11 




4 

0.0030 

cos 

(- g 





4 

«) 

12 




- 

0.0024 

cos 

(*2g 


♦ 

2F 


4 

ft) 

13 




4 

0.0023 

cos 

( + g 


♦ 

2F 


4 

0) 

14 




4 

0.0022 

cos 

(“ g 


♦ 

2F 

- 2D 

4 

2ft) 

15 




4 

0.0014 

cos 

« 

( 


4 

2F 

♦ 20 

4 

2ft) 

16 




' 

0.0011 

cos 

(♦ s 


4 

2F 

• 2D 

4 

2ft) 

17 




4 

0.0011 

cos 

( + 2g 


♦ 

2F 


4 

2ft) 

18 




•> 

o;ooio 

cos 

(- g 


4 

2F 


4 

a) 

19 




4 

0.0008 

cos 

( 


g 1 



4 

ft) 

20 




- 

0.0007 

cos 

(- g 

fcaua 



♦ D 

4 

ft) 

21 
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TABLE lb (Cont.) 





Series 

No. 

- 0.0007 cos 

(- g 

- 

2D + 

ft) 

22 

+ 0.0007 cos 

(+ g ♦ 

2g ’ + 2F - 

2r + 

2D) 

23 

+ 0.0005 cos 

C 

g* 

+ 

ft) 

24 

+ 0.0005 cos 

(- g 

+ 2F + 

2D + 

ft) 

25 

- 0.0003 cos 

( + 

g' + 2F 

+ 

2D) 

26 

+ 0.0003 cos 

( 

g * + 2F 

+ 

2ft) 

27 

+ 0.0003 cos 

O g 

+ 2F + 

2D + 

2ft) 

28 

♦ 0.0003 cos 

( 


2D + 

ft) 

29 

♦ 0.0003 cos 

C-2g 

♦ 

2D + 

0) 

30 

+ 0.0003 cos 

( 

g' + 2F - 

2D + 

ft) 

31 

- 0.0003 cos 

(+ g 

+ 2F - 

2D + 

ft) 

32 

+ 0.0003 cos 

( 

- 

2D + 

ft) 

33 

+ 0.0003 cos 

( 

+ 2F + 

2D + 

ft) 

34 

- 0.0002 cos 

( + 2g 

♦ 2F - 

2D + 

2ft) 

35 

+ 0.0002 cos 

( 

2g' + 2F - 

2D + 

ft) 

36 

- 0.0Q02 cos 

(*2g 

- 

2D + 

ft) 

37 i 

•t 0.0002 cos 

( + 2g 

+ 2F 

+ 

ft) 

38 

- 0.0002 cos 

c 

g’ ♦ 2F - 

2D + 

ft) 

39 

+ 0.0002 cos 

C-2g 

♦ 2F 

+ 

2ft) 

40 


TABLE lc NUTATION IN LONGITUDE 


Series No, 


♦ (-0V01737 T - 17V2327) sin ( 

♦ (-0.00013 T - 1.2729) sin ( 

+ (+0.00002 T + 0.200G) sin ( * 

♦ (-0.00002 T - 0.2037) sin ( 

♦ (-0.00031 T + 0.1261) sin ( 

+ (+0.00001 T + 0.067S) sin (♦ g 

+ (+0.00012 1 - 0.0497) sin ( 

+ (-0.00004 T •• 0.0342) sin ( 


+ 2F - 


♦ 2F 


g« + 2F 
♦ 2F 


♦ ft) 

2D ♦ 2 ft) 

♦ 2D) 

♦ 2fl) 


2D + 2ft) 

+ 2ft) 


13 





TABLE lc CCont.) 






Series 

No . I 


0.0261 

sin (+ g 


ft) 

9 

+ r-C, 00005 T + 

0.0214) 

sin ( 


g' + 2F - 2D + 2Q) 

10 

- 

0.0149 

sin (+ g 


- 2D ) 

11 

+ (+0.00001 T + 

0.0124) 

sin ( 


+ 2F - 2D + ft) 

12 

+ 

0.0114 

sin (- g 


+ 2F + 2D) 

13 

+ 

0.0060 

sin ( 


+ 2D ) 

14 

+ 

0.0058 

sin (+ g 


' + ft) 

15 

- 

0.0057 

sin (- g 


♦ ft) 

16 

- 

0.0052 

sin (- g 


+ 2F + 2D ♦ 2ft) 

17 


0.0045 

sin (-2g 


+ 2F + ft) 

18 

+ 

0.C04S 

sin (+2g 


- 2D ) 

19 

- 

0.0044 

siri (+ g 


+ 2F ♦ ft) 

20 

- 

0 .0032 

sin ( 


+ 2F + 2D ♦ 2ft) 

21 

+ 

0.0028 

sin (+2g 


) 

22 

+ 

0.0026 

sin (+ g 


♦ 2F - 2D ♦ 2ft) 

23 

- 

0.0026 

sin ( + 2g 


+ 2F + 2ft) 

24 


0.0025 

sin ( 


♦ 2F ) 

25 


.0.0621 

sin ( 


+ 2F - 2D ) 

26 

+ 

0.0019 

sin (- g 


♦ 2F ♦ ft) 

27 

♦ (-o/ooooi T + 

0.0016) 

sin ( 

♦ 

2g * ) 

28 

+ (+0.00001 T - 

0.0015) 

sin ( 


2g' ♦ 2F - 2D + 2ft) 

29 

- 

0.0015 

sin ( 

♦ 

g* ♦ ft) 

30 

+ 

0.0014 

sin (- g 


♦ 2D ♦ ft) 

31 

- 

0.0013 

sin (+ g 


- 2D ♦ ft) 

32 

- 

0.0010 

sin ( • 

- 

g' * 

33 

+ 

0.0010 

sin (+2g 


- 2F ) 

34 

- 

0.0009 

sin (- g 


2F + 2D + ft) 

35 

■f 

0.0007 

sin ( 

♦ 

g* * 2F •> 2ft) 

36 

- 

0.0007 

sin (+ g 

♦ 

ft - 2D ) 

37 

♦ 

0.0006 

sin (+ g 


+ 2D ) 

JL- 
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0.0006 

sin 

c 

g« + 

2F 


+ 

2D) 

39 


0-0006 

sin 

(+ g 

+ 

2F 

4* 

2D + 

20) 

40 

+ 

0.0006 

sin 

(+2g 


2F 

- 

2D + 

20) 

41 

- 

0 .0006 

sin 

c 



+ 

2D ^ 

0) 

42 

«* 

0.0005 

sin 

t*2g 




2D + 

0) 

43 

• 

0.0005 

sin 

c 

g’ + 

2F 

- 

2D + 

0) 

44 

4' 

0.0005 

sin 

( g 

♦ 

2F 

- 

2D + 

0) 

45 

- 

0.0005 

sin 

( 



m 

2D + 

0) 

46 

_ 

0.0005 

sin 

( 

♦ 

2F 

+ 

2D + 

0) 

47 

- 

0.0004 

sin 

( 

2g’ ♦ 

2F 

- 

2D + 

0) 

48 ] 

J 
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TABLE 2b NUTATION In OBLIQUITY 






Series 

No . 

A- - + 9V2100 

cos 

f 


+ ft) 

1 

+ 0.5522 

cos 

( 


+ 2F - 2D + 2ft) 

2 

- 0.0904 

cos 

( 


* 2ft) 

3 

+ 0.0884 

cos 

( 


+ 2F + 2ft) 

4 

+ 0.0216 

cos 

( 

g’ 

+ 2F - 2D + 2ft) 

5 

+ 0.0185 

cos 

r 


♦ 2F ♦ ft) 

6 

+ 0.0113 

cos 

( + g 


+ 2F + 2ft) 

7 ' 

- 0.0093 

cos 

( _ - 

g* 

+ 2F - 2D + 2ft) 

8 

- 0.0066 

cos 

( 


+ 2F - 2D + ft) 

9 

- 0.0050 

cos 

C- g 


♦ 2F + 2ft) 

10 

- 0.0031 

cos 

( + g 


♦ ft) 

11 

+ 0.0030 

cos 

(- g 


+ ft) 

12 

- 0.0024 

cos 

(-?g 


+ 2F + ft) 

13 

+ 0.0023 

cos 

( + g 


+ 2F + ft) 

14 

+ 0.0022 

cos 

(" g 


+ 2F + 2D + 2ft) 

15 

+ 0.0014 

cos 

( 


+ 2F + 2D + 2ft) 

16 

- 0.0011 

cos 

( + g 


+ 2F - 2D + 2ft) 

17 

+ 0.0011 

cos 

( + 2g 


+ 2F + 2ft) 

18 

- 0.0010 

cos 

(* g 


+ 2F + ft) 

19 




TABLE 2c NUTATION IN LONGITUDE 


* + (-O'. '01737 T - 



] 772327) siii 
1.2729) sin 
0.20QS) sin 


+ 2F 


Series No. 

+ ft) 1 

2D + 2D) 2 

+ 2D) 3 


0.2037) 

s j n 

( 


+ 

2F 


4 

2D) 

4 

0.3 261) 

sin 

c. 

g* 





) 

5 

0.0675) 

sin 

( 4 e 






) 

<r 

Sj 

0.0497) 

sin 

( 

g* 

4 

2F - 

2D 

+ 

2D) 

7 

0.0342) 

sin 

c 


4 

2F 


4 

2D) 

8 

0.0261 

sin 

( + R 


+ 

2F 



a) 

9 

0.0214) 

sin 

( 

c' 

4 

2F - 

2D 

V 

2D) 

10 

0.0149 

sin 

C 4 g 



- 

2D 


) 

11 

0.0124) 

sin 

c 


4 

2F - 

2D 

- * 

D) 

12 

0.0114 

sin 

C- g 


4 

2F 


4 

2D) 

13 

0.0060 

sin 

c 



+ 

2D 


) 

14 

0.0058 

sin 

C + g 





4 

ft) 

15 

0.0057 

sin 

(- g 





4 

ft) 

16 

0.0052 

sin 

(*• g 


4 

2F + 

2D 

4 

2ft) 

17 

0.0045 

sin 

c- 2 g 


4 

2F 


4 

ft) 

18 

0.0045 

sin 

( + 2g 



- 

2D 


) 

19 

0.0044 

sin 

( + g 


4 

2F 


4 

o) 

20 

0.0032 

sin 

c 


4 

2F 

2D 

4 

2ft) 

21 

0.0028 

sin 

("•2g 






) 

ll 

0.0026 

sin 

( + g 


4 

2F - 

2D 

4 

2ft) 

23 

0.0026 

sin 

(*2g 


4 

2F 


4 

2ft) 

24 

0.0025 

sin 

c 


4 

2F 



) 

25 

0.0021 

sin 

( 


4 

2F - 

2D 


) 

26 

0,0019 

sin 

(- g 


■ 4 

2F 


4 

ft) 

27 

0.0016) 

sin 

( 

2g 

♦ 




5 

28 



TABLE 2c (Cont.J 


0.0015; 

sin 

c 


2g* + 2F - 2D 

4 

2D) 

29 

C . 0015 

sin 

c 

+ 

g’ 

+ 

to 

30 

o 

o 

o 

sin 

c- 

g 

♦ 2D 

4 

fl) 

31 

0.0013 

sin 

o 

g 

- 2D 

4 - 

ft) 

32 

0.0010 

sin 

( 

- 

g* 

4 

to 

33 

0.0010 

sin 

( + 

2g 

- 2F 


) 

34 


To eh* time units for coefficient of 1st term, use 
.475 S6i> '*0* 6 d - .01757 T 


Table 2d: True obliquity of Date and Nutation in -'ight ascension 


c T * e M ♦ 

A a * cos 


t 


i 

ii 
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SECTION 4.0 

LUN I - SOLAR - PLANET.' RY EPHEMERIS 

GEODYN uses precomputed equi-speccd ephemeris data EPHEM 
in true of date coordinates for the Moon, the Sun, Venus, 

Mars, Jupiter ar.d Saturn. The actual ephemerides are com- 
puted using Everett's fifth-order interpolation formula. 

The interval between ephemerides; i.e., the tabular interval 
h, is 0.5 days for the Moon and the Equation of the equinoxes 
and 4.0 days for the other bodies. 

The GEODYN ephemeris tape contains pH coordinates 
in true of date. The quantities on the tape are 

a) geocentric lunar positions and the corresponding 
2nd and 4th dif lure aces , 

b) solar positions relative to the earth-moon 
barycenter and the corresponding 2nd and 4th 
differences , 

c) heliocentric positions of Venus, Mar c , Jupiter 
and Saturn and the corresponding 2nd and 4th 
differences , 

d) - the equation of the equinoxes and its 2nd and 

4th differences. 

The format of this tape is presented in Volume III 
of the GEODYN documentation 

This ephemeris tape was prepared from a JPL planetary 
ephemeris tape corresponding to "JPL Development Ephemeris 
Number 69," Reference 1. The program which generates the 
GEODYN ephemeris tape is described in Volume IV of the 
GEODYN documentation. 


. 0-1 


y (t ^ + sh) 


F 0 (l-s)*d j 2 F 2 (1 -s) 


* i. A F 4 (1-s) 


♦*,♦1 V s ^ d jn V s) 


•» d j+1 F 4 (s) 


F„'S) = S 


F 2 (s) » 1(S-1) (S) (s+11]/6 

F 4 (s) = [ (s - 2) (s-1) (s) (s+1) (s+2j]/120 


The quantity s is of course the fractional interval for 
the interpolation. The quantities are obtained from 
the ephemeris tape. 


t i 


The formulation for Everet-’s fifth-orde- interpo- 


lation is 


where 


EPHEM 
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SECTION 5.0 
THE OBSERVER 


This section is concerned with the position and 
coordinate systems of the observer. Thus it will cover 

• geodetic station position coordinates, 

• topocentric coordinate systems, 

• time reference systems, and 

• polar motion. 

The geodetic station position coordinates are a 
convenient and quite common way of describing station 
positions. Consequently, GEODYN contains provisions for 
converting to and from these coordinates, including the 
transformation of the covariance matrix for the deter- 
mined Cartesian station positions. 

The topocentric coordinate systems are coordinate 
systems to which the observer references his observations. 

The time reference systems are the time systems in 
which ihe observer specifies his observations. The 
transformations between time reference systems are also 
given. These latter are used both to convert the ob- 
servation times to A1 time, which is the independent 
variable in the equations of motion, and to convert the 
GEODYN output to UTC time, which is the generally recognized 
system for output. 
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The positions of the observers in GEODYN are referred 
to an Earth-fixed system defined by the mean pole of 1900.5 
and Greenwich. They are rotated into the Earth-fixed system 
of date at each observation time by applying "polar motion", 
which is considered to be slippage of the Earth’s crust. 


5.1 GEODETIC COORDINATES 

Frequently, it is more convenient t- define the 
station positions in a spherical coordinate system. 

The spherical coordinate system uses an oblate spheroid 
or an ellipsoid of revolution as a model for the geo- 
metric shape of the Earth. The Earth is flattened 
slightly at the poles and bulges a little at the 
equator; thus, a cross section of the Earth is approxi- 
mately an ellipse. Rotating an ellipse about its 
shorter axis forms an cblate spheroid. 

An oblate spheroid is uniquely defined by specify- 
ing two dimensions, conventionally, the semi-major axis 
and the flattening, f, where f * (dee Figure 1) 

a 

This model is used in the GEODYN system. The 
spherical coordinates utilized are termed geodetic co- 
ordinates and are defined as follows: 

• $ is geodetic latitude, the acute angle 

between the semi -major axis and a line 
through the observer perpendicular to 
the spheroid. 


3 X is east longitude , the angle measured 

a rh& r»nn«tnrial Diane between 

w>a j i tju* *> ■> — * n ' * * 

the Greenwich meridian and the observer s 
meridian. 

• h is spheroid height, the perpendicular 

height of the observer above the refer- 
ence spheroid. 

Consider the problem of convert^* from X, and 

h to X , Y , and Z * the Earth- fixed Cartesian coordinates. 
© © ^ • 

The geometry for an X-Z plane is illustrated in 
Figure 1 . The equation for this ellipse is 

X 2 ♦ * 5- ■ a 2 » ^ 

( 1-0 

where the eccentricity has been determineo om the 
flattening by the familiar relationship 


e 2 - 1 - Cl-f) 2 . 


( 2 ) 


SQUANT 


5 . 1-2 




-+t> 


The equation for the normal to the suriace ui 
eiipse yields 


c r ri hvt 


tan <j> 


By taking differentials on equation (1) and applying 
the result in equation (3), we arrive at 


— 3 ( 1 -e ) tan 4 

X 


The simultaneous solution of equations (1) and (4) for 
X yields 


a cos <J> 


' 1 -e sin 


From inspection of Figure 1 we have: 


cos <j> * - ; 

N 


and hence, applying equation (5), 


'1-e sin 
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For an observer at a distance h from the refer- 


SQUANT 


ence ellipsoid, the observer’s coordinates (X,Z) become 


X » N cos <}> + h c.ts 4> (8) 

and 

Z * N (1-e 2 ) sin 4> + h sin (9) 

The conversion of <J>, X, and h to X . Y . and Z 

is then 


"V 


”(N+h) cos <{> cos X - 


s 

(N+h) cos 4> sin X 

. Z o_ 


(N+h-e 2 K) sin 4> 


In the GEQDYN system this conversion is performed in 
subroutine SQUANT. 

The problem of converting from X g , Y , and to 
4, X, and h is more complex as we cannot staTt with a 
point on the reference ellipsoid. For this reason the 
determination of accurate values for $ and h requires 
an iterative technique. 


Conversion to Geodetic Coordinates 


For the problem of converting station coordinates 
in X , Y . and to 4>, X, and h we know that N is on the 
order of magnitude of an Earth radius, and h is a few 
meters. Hence 


h << N 


( 11 ) 


The Earth is approximately a sphere, hence 


e << 1. 


( 12 ) 


Therefore, again working in our X-Z plane (see Figure 1) , 


N sin $ ~ Z. 


(13) 


From Figure 1 (see also equation (9)) we have 


t 


Ne 2 sin $, 


(14) 


or, for an initial approximation, 


t 


Z. 


(IS) 
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PLHOUT 


PLHOUT 


( 


The series of calculations to be performed on 
each iteration is: 


Z t = z ♦ t 


N*h - ♦ y 7 ♦ i\ j 


, \ 1/2 


sin ■ Z / 


(N+h) 


( 


" * V(l-e J sin **) 1 ' 2 


t * Ne sin <j>. 


( 16 ) 

( 17 ) 

( 18 ) 

( 19 ) 

( 20 ) 


When t converges, 4> and h are computed from sin and 
(N+h) . The computation of \ is obvious; it being simply 


X ‘ * 



( 21 ) 


This procedure for determining X, and h is that coded 
in subroutine PLHOUT. 
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There is a different procedure in subroutine 
PREDCT for computing <j>, X, and h for a satellite. This 
is because the accuracy requirements are less stringent. 

This different procedure it also used in subroutine 
BRAG to evaluate the satellite height for subroutine DENSTY. 

Because e << 1, '-e may write an approximation 
to equation (9) : 


1 « (N+h) (1-e 2 ) sin 4> = ' (22) 

From Figure 1, 


X 


(N+h) cos <}> 



+ 



(23) 


and by remembering equation (2), 



tan 


Z 


e 




(24) 
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PREDCT 


The equation for the ellipse, equation Cl) , 
yields the following formula for the radius of the 
ellipsoid: 


ellipsoid 


P~T?. 


a (1-f) 


/l'- C2f-f 2 ) "(l-sin^T") 


(27) 


where <j>" is the geocentric latitude. After applying the 
Binomial Theorem, we arrive at 


ellipsoid 


[ 


a I 1 - (f+| f 2 ) sin 4 <j>' + £ f 4 sin* 


3 <-2 ..4 


C28) 


wherein terms on the order of f have been neglected. The 
(spheroid) height may then be calculated from r, the geo- 
centric radius of the satellite: 


r - r 


ellipsoid, or 


(29) 


V s 


~ 1 — 7 — T 


h . - */Xg+ Yg+Zg - a + (af+|af 2 ) sinV -|af 2 sinV (30) 


The sin| of the geocentric latitude, sin is of 
course -=r-. 

T 
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Subroutine VEVAL also requires the partial 
derivatives of h with respect to position for the diag 
variational oartials computations: 


VEVAL 


+ 2 sin 


I / 3 

s v [af + ~ af* 


3 af 2 sin 2 $ 


2 r. 13: 


e 1 


r 3r . 


(31) 


where the 


are the Earth- fixed components of r; i.e., 

V Z e>- 


In addition to the conversion of the coordinates 
themselves, GEODYN also converts covariance matrices for 
the station positions to either the 4», X, h system or 
the Earth- fixed rectangular system. This is accomplished 
in INOUPT, SQUANT, and PLHOUT by calling VCQNV to compute 


INOUPT 

SQUANT 

PLHOUT 

VCONV 


V OUT ' PTV IN P 


1321 VCONV 


where V Q y T is the output covariance matrix, Vj N is the 
input covariance matrix, and P is the matrix of partials 
relating the coordinates in the output system to the 
coordinates in the input system. 


I! I 
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These partial derivatives (ir P) which GEODYN PLHOUT 

requires are for Y , with respect to <J> . X, h and 
vice versa. These partials are: 


-X e Z e (l-e 2 )/((l-e 2 ) 2 (X^Yj) + Z 2 ) (x 2 +y < 2,1 


3X C 


e e-' e / v ‘ e e 


§4 - 'Y e Z e (l-e^/C(l-e 2 ) 2 (X 2 + Y 2 ) + Z 2 ) (xf + Y 2 ) 2 


3X 


|| - (X 2+ Y 2 ) (l-e 2 )/(l-e 2 ) 2 (X 2 +Y 2 ) + Z 2 ) (X 2 +Y 2 ) 2 


e e 


e 


e e J 


V!,‘ 

It’ 


(33) 


3X 

7T 


It = fx C-e 2 a.(l-e ? ') s iT\<h cos*/(l-e 2 sin 2 <j>)^ -Z e cos<j>/sin 2 <t>) 


3h 3^ 
3Y e 


(-e 2 a(I-e 2 )sin<}> cos<j>/ (1- e 2 sin 2 4>)^ -Z e cos<|>/sin 2 <{>) 


|| - || (-e 2 a(l-e 2 )sin<t> cos<J>/ (1-e 2 sin 2 <j>) ? * Z e cos<t>/sin 2 ♦) 
& e 


sa.ii <{> 
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f 


f 


I 


1 


I 

I 



-sin$ cosA \ N+h 

L 


2 2 
No cos ‘b 

7 *“ - 

1-e sin 


a *e 

-(N+h) co'*'}' sinX 

3X 



cos4» cosX 



r 


-sin<{> sin\ 


N+h 


2 2 
Ne cos £ 4> 

' T . 2 7 

1-e sin <|> 



f 

i 


(N+h) c^-<)) cosX 


(34) 


( 


3h 


* cos<{> sinX 


aZ. 




= cos<J> 


e^sin^ <k 

h+N (1-e") \ 1+ 5 T ~ 

' 1-e siii <? 


az. 


ax 


5 

3h 


= sin<j> 



The partials for converting from X .Y . to 

v C C 

X, h are computed in subroutine PLHOUT. Those for 
converting fro»n <J>* X, h to X ,Y ,Z are computed in 
subroutine SQUANT. 


SQUANT 
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To:-v::L.aRic coordinate systems 






i 






( 



The observations of a spacecraft are usually 
ief^i-uccJ to the observe., and therefore an additional 
set of reference systems is used for this purpose. The 
origin of these systems, referreu to as topocentric 
coordinaLC systems, is the obsc-'"’ :i on the surface of 
•'■he earth. 


Topocentric right ascension and declination are 
measured in an inertial system whose Z axis and funda- 
mental plane are parallel to those of the geocentric 
inertial *>stem. The X axis in this case also points 
toward the vernal equinox. 


Th other major topocentric system is the Earth- 
system determined by the zenith and the observer's 
horizon plane. This is an orthonormal system defined 

A A A 

by N, E, and Z, which are unit vectors which point in 
the same directions as vectors from the observer 
pointing north, east, and toward the zenith. Their 
definitions are: 


N 


- sin 4 cos 

- sin $ sin X 

cos $ 


E - 


sin ' 
ccs X 
o 


A 

7 

u 


cos 4> cos X 
cos 4* sin X 
sin 4 


( 1 ) 


( 2 ) 


( 3 ) 
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SQUANT 


SQl'ANT 

PREDCT 

OBSDOT 


where <fc is the geodetic latitude and X is the east 
longitude of the observer (see Section S.l). 

A A A 

These N, E, and Z vectors are computed in 
SQUANT for use in PREDCT and OBSDOT. 

This latter system is the one to which such 
•measurements as acinuth and elevation, X and Y angles, 
and direction cosines are related. 

It should be noted that the reference systems for 
range and range rate must be Earth- fixed, but the choice of 
origin is arbitrary. In GECDYN , range and range rate are 
not considered to be topocer.tric, but rather geocentric. 


3.3 TIME REFERENCE SYSTEMS 

Three principal time systems are currsn..; in 
u ephemeris time, atomic time, and universal time. 


Ephemeris time is the independent variable in 
the equations of motion of * e sun; this time is the 
uniform mathematical time. The corrections that must 
be applied to universal time to obtain ephemeris time 
are published in the American Ephemeris a>i Nautical 
Almanac or alternatively by BIH, the "Bureau Inter- 
national de l'deure." 

Atomic time is a time based on the oscillations 
of cesium at zero field. In practice A1 time is based 
on the mean frequency of oscillation of several cesium 
standards as compared with the frequency of ephemeris 
tire. This is the time system in which the satellite 

equations of motion are integrated in GEODYN. 
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Universal time is determined by the rotation of 
'•’•e Earth. UT1 . the time reference system used ir* 
GEODYN to position the Earth, is universal time that 
has been corrected for polar motion. UTC is the time 
of the transmitting clock of any of the synchronized 
transmitting time signals . The frequency of a UTC 
clock is pre-set to a predicted frequency of UT2 time, 
where UT2 time is universal time corrected for ob- 
served polar motion and extrapolated seasonal variation 
in the speed of the earth’s rotation. 

The reader who is unfamiliar with these time 
systems should refer to one of the annual reports o£ 
BIH. 


5.3.1 Time System Transformations 

The time system transformations are between any 
combination of the A1 , UT1, UT2, or UTC referv sys- 
tems. These transformations are computed in the 
GEODYN system by subroutine TDIF. 

The time transformation between any -nput time 
system and any output time system is formed by simple 
addition and subtraction of the following set of time 
differences : 

• UT2 - UT1 

• A1 - UT1 

e A1 - UTC 
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TDiF 


TDIF 


The following equation is used to calculate 
(UT2-UT1) for any year: 

(UT2-UT1) = + *022 sin 2*,tt-!oi2 cos 2irt CD 

- ^006 sin 4Trt+!o07 co c ... c 

t » fraction of the tropical year 

elapsed from the beginning of the 
Besseliau year for which the 
calculation is made. 

(1 tropical year = 365.2422 days) 

This difference, (UT2-UT1) , is also known by the name 
'’seasonal variation." 

The time difference (Ai-UTl) is computed by 
linear interpolation from a table values. 

The spacing for the table is every 10 days, which 
matches the increment for the "final time of emission" 
data published by the U.S. Naval Observatory in the 
bulletin, "Time Signals," The differences for this 
table are determined by 


(A1 - ’JT1) » (A1 - UTC) - (UT1 - UTC) 

The values for (UT1 - UTC) are obtained from "Circula- D", 
B1H. The differences (A1 - UTC) are determined according 
# to the following procedure. 
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The computation of (Al-UTC) is simple, but net 
so straightforward, UTC contains discontinuities both 
in epoch and in frequency because an attempt is made 
to keep the difference between a UTC flock and a UT2 
clock less than .1. When adjustments are made, by 
international agreement they are made in steps of ?1 
and only at the beginning of the month; i,e., at o*?c UT 
of the first day of the month. The general formula which 
is used to compute (Al-UTC) is 


(Al-UTC) = a Q + a x (t-t Q ) 


( 2 ) 


Both a,p and a n are recovered from tables. The values 
in the table for a Q are the values of [Al-UTC) at the 
time of each particular step adjustment. The values 
in the table for a^ are the value? for the new rates 
of change between the two systems after each step 
adjustment . 

Values for a Q and a, are published both by the 
U.S. Naval Observatory and BIH. 


S.4 POLAR MOTION 

Consider the point P which is defined by the 
intersection of the Earth's axis of rotation at some 
time t with the surface of the Earth. At some time t+At, 
the intersection will be at some point P' which is diff 
than P. Thus the axis of rotation appears to be moving 
tiv? to a fixed position on the Earth; hence the term " 
of the pole." 

5.4-1 

• - ■ » i 



POLE 


POLE 


Le us establish a rectangular coordinate sys 1- 

O 

«*(»ntered at a point F fixed u n the surface of the Ea. a 
with F near the point P around 1900, and take measure- 
ments of the rectangular coordinates of the point P 
during the period 1900.0 - 1906.0. It is observed that 
the point P moves in roughly circular motion in this 
coordinate system with two distinct periods, one period 
of approximately 12 months and one period of 14 months. 

We define the mean position of P during this period to 
be the point Pq , the mean pole of 1900.0 - 1906.0. 

The average is taken over a six year period in 
order to average out both the 12 month period and the 
14 month period simultaneously (since 6 times 12 months * 
72 months and 72/14 = 5 periods approximately of the 
14 month term) . The radius of this observed circle 
varies between 15-35 feet. 

In addition to the periodic motion of P about P Q , 

by taking six year means of P in the years after 1900 - 

1906, called P , there is seen to be a secular motion 
m 

of the mean position of the pole away from its original 
mean position Pq in the years 1900 - 1906 at the rate of 


» 
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approximately 0'.'0032/year in the direction cf the POLE 

mcridan 60° W, and a libration motion of a period of 
approximately 24 years with a coefficient of about 
07022. The short periodic motions over a period of 
six years average about 072 - 073. 


5.4.1 Effect on the Position of a Station 


This motion of the pole means that the observing 
stations are moving with respect to our "Earth- fixed" 
coordinate system used in GEODYN. The station positions 
must be corrected for this effect. 

The position of the instantaneous or true pole 
is computed by linear interpolation in a table of ob- 
served values for the true pole r'lative to the mean 
pole of 1900 - 1905. The table increment is 10 days; 
the current range of data is from December 1, 1960 to 
June 1, 1972. The user. should be aware of the fact 
that this table is expanded as n*.w information becomes 
available. If the requested time is hot in the range of 
the table, the value for the closest time is used. 

The data in the table is in the form of the co- 
ordinates o^ the f*ue pole relative to the mean pole 
measured in seconds of arc. This data was obtained from 
"Circular D" which is published by BIH. The appropriate 
coordinate system and rotation are illustrated in Figures 
1 and 2. 





st uLbas* ax " 




I 



= Center of Coordinate System 
= Adopted Mean Pole 

X-j = Direction of 1 st Principal Axis (along meridian 
directed to Greenwich) 

X 2 = Direction of 2 n< * Principal Axis (along 90° 

West meridian) 

Py * Instantaneous Axis of Rotation - 
x,y * Coordinates of ?y Relative to P^ Measured 
in seconds of arc 

Figure 1: Coordinates of the Instantaneous Axis of Rotation 


*3 

V Y 3 y x p a 


X.y^r 


A^.y 2 


4 

x,y * Rectangular Coordinates of Py Relative to P^ 

X,X 2 Plane « Mean Adopted Equator Defined by 
Direction of Adopted Pole P A 

Ya/ 2 Plane * Instantaneous Equator Defined by 
y c Direction of Instantaneous Pole Py 

Rotation of Coordinate System from Adopted Mean Pole 
System to Instantaneous Pole System 
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Consider the station vector X in a system attached 
to the Earth of the mean po"! e and the same vectoi V 
in the "Earth- fixed" system of GEODYN . The transforma- 
tion between Y and X consists of a rotation of x about 
the X 2 axis and a rotation of y about the X^ axis; 
that is 


R-j (y) R 2 (x) X 


0 cos y sin y 

0 -sin y cos y 


cos x 


sxn x 


sin x 


cos x 


Because x and y are small angles, their cosines 
are set to 1 and their sines equal to their values in 
radians. Consequently, 


In the GEODYN system, the position of the true 
pole is computed by subroutine POLE. The station vec- 
tors are referenced to the true pole by subroutine 
TRUEP. 


POLE 

TRUEP 
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The coordinate rotation is defined as 


TRU12P 


u * -^(y) w 


PRHDCT 

CD .. 


where 

* ♦. 

w -• station vector in a system attached to the 
Earth of the me„ ■* pole, 

# 

% 

u = station vector in a system attached to the 
Earth of the true pole. 

R 1 (y) = matrix of rotation about the X^ axis 

« * 

R^,(x) = matrix Qf rotation about the X 2 axis 
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The rotation matrices are. 


R x iy) 


10 3 

0 cos y sin y 

j -sin y cos y 

- 


R 2 U) 


cos x 0 

0 1 

sin x 0 


-sin x 

0 

cos x 


Defining 


u = u, i + u ? j + u, k 


17 = w, i + w~ j + w, k 


( 2 ) 

( 3 ) 


and performing the matrix multip.’ icatiqns . 


u, « w, cos x - sm x 
1 1 i 


u 7 * s? n x sm y + cos y + cos sin y 


( 4 ) 


u 3 * "1 


w, sir x cos y - w~ sin y + w_ cos x cos y 
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The fundamental quantities required for the esti- 
mation of polar motion parameters are 


(fill 0m. 

— • and — 

ax by 


where m is the satellite observation, 
Using the chain rule 


PREDCT 


dm 

3m 

3u^ 3m 

3u 2 3m 

SUj 



3x 3«2 

ax auj 

3x 

3m 

dm 

3Uj 3in 

9 ii 2 am 

9u 3 

ay 

3u i 

3 y 3u^ 

3y 3v 
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The partial clcrvia t. ivcs oi the satellite observation 
with respect to the true station coordinates are currently 
available in CEODYN, The partial derivatives of the 
station coordinates with respect to the polar motion 
parameters are: 


- -w^ sin x 


Wj cos x 


• • - 0 


— — = w. cos x sin y - w- sm x sir* y 

3x 1 ° 


= w. sin x cos y - v ?2 s 4 , y + w^ cos x cos y 


= w^ cos x ccs y 


Wj sin x cos y 


= -w^ sin x sin y - Wj cos y - w^ cos x sin y 


? . 4 - 9 


r 


"T 
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Since the angles x and y are small, the following 
approximations may be made. 


sin x = x 


cos x * 1 


sm y - y 


cos y - 1 


Suds ituting equations (7) into equation*- (6) 


TRUEP | ! 


w„ x - w, 

1 3 


= w, y - w, x y 


w i x “ '2 y + w 3 


* W 1 “ w 3 ** 


‘ w i x y ' w 2' ~ w 3 y 
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MEASUP -MENT MODELING ANT RELATED DERIVATIVES 


The observations in GEODYN are geocentric in nature 
The computed values for the observations are obtained by 
applying these geometric relationships to the computed 
values for the relative positions and velocities of the 
satellite and the observer at the desired time. 

In addition to the geometric relationships GEODYN 
allows for a timing bias and for a constant bias to be 
associated with a measurement type from a given station. 
Both of these biases are optional. 

The measurement model for GEODYN is therefore 
• • * 

c t*4t " f t (7 - F > F ob> * b * f t (F - F « F ob )-it (1) 

where 

^t+At is the computed equivalent of the ob- 
servation taken at time t+At, 

r is the Earth-fixed position vector of 

the satellite. 

\ 9 

* 

is the Earth-fixed position vector of 
the station, 



f t (.r ,r ,r o ^j is liie geometric relationship defined 

by the particular observation type at 
time t, 

b is a constant bias on the measurement, 

and 

At is the timing bias associated with the 

measurement . 

The functional dependence of f^ was explicitly stated for 
the general case. Many of the measurements are functions 
onl) of the position vectors and are hence not functions 
of ;he satellite velocity vector r. We will hereafter :*cfer 
to f t without the explicit functional dependence for rota- 
tional convenience. 

As was indicated earlier in Section 2.2, we require 
the partial derivatives of the computed values for the 
measurements with respect to the parameters being determined 
(see also Section 10.1). These parameters are: 

• the -rue of date position and velocity of the 
satellite at epoch. These correspond to the 
inertial position and velocity which are the 
initial conditions for the equations of motion, 

• force model p«nameters, 

• the E:rth-f ixed station positions, 


• measurement biases. 


Thes>? parameters are implicitly divided into a set 
a which are not concerned with the dynamics of satellite 
motion, and a set F which are. 

The partial derivatives associated with the param- 
eters a; i.e., station positions and measurement biases 
are computed directly at the given observation times. The 
partial derivatives with respect to the parameters F; i.e., 
the epoch position and velocity and the force model param- 
eters, roust be determined according to a chain rule: 


3C t+At = 3C t+At 
3F 3x t 

where 

x t is the vector which describes the satellite 
position and velocity in true of date co- 
ordinates. 

3C t+At 

The partial derivatives — — are computed directly at the 

given observation times, but the partial derivatives — - 

3F 

may not be so obtained. These latter relate the true of 
date position and velocity of the satellite at the given 
time to the parameters at epoch through the satellite 
dynamics . 


3x. 


33 


( 2 ) 
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The partial derivatives are called the varia- 

aF 

tional partials and are obtained by direct numerical 
integration of the variational equations. As will be 
shown in Section 8.2, these equations are analogous 
to the equations of motion. 

Let us first consider the partial derivatives of 
the computed values associated with the parameters in 6. 
We have 


SC t+At m Ht ^t ( 

3F 3x t 30 

Note that we have dropped the partial derivative with 

* • 

respect to F of the differential product f^At. This is 
because we use first order Taylor series approximation 
in our error model and hence higher order terms are 
assumed negligible. This linearization is also com- 
pletely consistent with the linearization assumptions 
made in the solution to the estimation equations 
(Section iO. 1) » 

3f t 

The partial derivatives are computed by 

3x t 3£ t 3f t 

transforming the partial derivatives and — *-=■ 

3r 3r 

from the Earth- fixed system to the true of date system 
(see Section 3.4). These last are the partial deriva- 
tives of the geometric relationships given later in this 
section (6. 2) . 


In summary, the partial derivative* icquired for 
3C . 

computing the -== , the partial derivatives of the 

3 6 

computed value for a given measurement, are the variational 
partials and the Earth-fixed geometric partial derivatives. 

The partial derivatives of the computed values with 
respect to the station positions are simply related to 
the partial derivatives with respect to the satellite 
position at time t: 


3C 


t+At 





( 4 ) 


where r is of course the satellite position vector in 
Earth- fixed coordinates. This simple relationship is a 
direct result of the symmetry in position coordinates. 
The function f is a geometric function of the relative 
position; i.e., the differences in position coordinates 
which will be the same in any coordinate system. 


The partial derivatives with respect to the biases 
are obvious: 


3C 

3b 


t+At 


(S) 


* C t+At 
3 (At) 


( 6 ) 
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In the remainder of this section, wc will be con- 
cerned with the calculation of the geometric function 
and its derivatives. These derivatives have been 
shown above to be the partial derivatives with respect 
to satellite position and velocity at time t and the 
time rate of change of the function, f t « 

The subroutine breakdown for the calculation of PREDCT 

these quantities in GEODYN is as follows: The geometric OBSDOT 

relationships and the geometric partial derivatives aTe 
implemented in subroutine PREDCT. The time rates of 
change are coded in subroutine OBSDOT. 

The data preprocessing also requires some use PROCES 

of these formulas for computing measurement equiva- 
lents. These are then also implemented in subroutine 
C PROCES. 



6.1 


THE GEOMETRIC RELATIONSHIPS 


The basic types of observation in GEODYN are: 

• riy.ht ascension and declination 

• range 

• range rate 

• l and m direction cosines 

• X and Y angles 

• azimuth and elevation 

• altimeter height and rate 

The geometric relationship which corresponds to each of 
these observations is presented below. It should be noted 
that in addition to the Earth- fixed or inertial coordinate 
systems, some of these utilize topocentric coordinate 
systems. These last are presented in Section S.2. 





Range : 


Consider the station-satellite vector: 


P 


* r - 



where 


( 1 ) 


r is the satellite position vector ( x,y,z ) in 
the geocentric Earth-fixed system, and 

r Q b is the station vector in the same system. 

The magnitude of this vector, p, is the (slant) 
range, whi'* u is one of the measurements. 

( 

Range rate: 

The time rate of change of this vector p" is 

i = 7 (2) 

1 as the velocity of the observer in the Earth-fixed sys- 

tem is zero. Let us consider that 

{) ■ pu (5) 

where 

u is the unit vector in the direction of 

» t 
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f 

GRHRAN 

PREDCT 

OBSDOT 

f 


Thus we have 


GRHRAN 

PREDCT 

OBSDOT 


• • A A 

F “ pu + pu 


( 4 ) 


The quantity p in the above equation is the computed value 
for the range rate and is determined by 



(5) 


Altimeter height: 

The altimeter height and rate are unique in that the 
satellite is making the observation. While these are 
actually measurements from the satellite to the surface 
of the Earth, they are taken to be measurements of the 
spheroid height and the time rate of change of that 
quantity for obvious reasons. Using the formula for 
spheroid height previously determined in Section 5.1. 
we have: r ( 



( 6 ) 


PREDCT 


6.1-5 


where 


PREDCT 


V5 


is the earth's mean equatorial radius, 

is the Earth's flattening, and 

is Tj, the z component of the Earth-fixed 
satellite vector. 


Altimeter rate: 


The altimeter rate is determined by a chain rule: 


PREDCT 


™alt ' 7 


The required partial derivatives are given in the section 
on geometric, partials. 
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Right ascension and declination : 


The topocentric right ascension a and declination 
6 are inertial coordinate system measurements as illus- 
trated in Fi/ure 1. GEODYN computes these angles from 
the components of the Earth-fixed station-satellite vec- 
tor and the Greenwich hour angle £>g. 



The remaining measurements are in the topocentric 
horizon coordinate system. These all require the N, Z, 

A 

and E (north, zenith, and east base line) unit vectors 
which describe the coordinate system. 




n; ~ ^ + 


recticn cosines 


There are three direction cosines associated with 
the station-satelli e vecto v in the topocentric system. 
These are: 


PREDCT 


— u • E 


‘The l and m direction cosines are observation types for 
GEODYN. 


X and Y angles: 

The X and Y angles are illustrated in Figure 2 
They are computed by 


3^ ■ tan 


-C-) 


sin’* (m) 


6.1-7 


I 

r 


Tracking Station 


z 

4 — . Lotol Horizontci Ptone 


FIGURE 2. X and Y Angles 




Figure 7 .illustrates the measurements of azimuth 
and elevation. These angles are computed by: 
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A = tan 
2 


E 4 = sin' 1 (n) 


(14) 


6.2 the GEONETRIC PARTIAL DERIVATIVES 


PREDCT 


The partial derivatives for each of the calculated 
geometric equivalents with respect to the satellite positions 
and velocity are given here. All are in the geocentric, 
Earth-fixed system. (The r^ refer to the Earth-fixed 
components of r . ) 

Range : 


Range rate: 
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Altimeter range; 
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X and Y Angles: PREDCT 
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3r i p(l-m^) 
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N^-rau^ 

p/l-m^ 


(17) 


Azimuth and Elevation: 
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mE.-tN. 

l i 

(18) 
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t 
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p/l-n^ ' 
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The derivatives of each measurement type with 
respect to time is presented below. All are in the 
Earth-fixed system. 

Range : 


u • r 


( 1 ) 


Range Rate : 

The range rate derivative deserves special atten- 
tion. ReveJnbering that 


( 2 ) 


We write 



u • p 


( 3 ) 


Thus 


u . p ♦ u • p 


( 4 ) 


.9 

„ ? 


OBSDOT 


I 


Because 


OBSDOT 


p = — (pu) = pu + pu 

dt 


we may substitute in Equation 4 above for u: 


•• 1 I • • A • ^ H 

p« - (p • P - p u • p) + u-p 
p 


or, as 


p » u • p 


we may write 


•• 1 • > • A •• 

P " - (P ** P * P ♦ P ‘ P) 


In order to obtain p, we use the limited gravity potential 
(see Section 8.3). 


« L . 

r \ 


~ ( sin ♦> 
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The gradient of this, potential with respect to the Earth- OBSDOT 
fixed position coordinates of the satellite is the part of 
"p due to the geopotentia! : 


3 a e C 20 


5 sin 4 <()- 1-2 


We must add to this the effect of the rotation of the 
coordinate syster . (The Earth- fixed coordinate system 
rotates with respect to the true, of date coordinates with 
a rate 0 g , the time rate of change of the Greenwich hour 
angle.) 


The components of p are then 


+ [x cos 6 g + y sin 0 g ] 6 g ♦ r 2 0 (11) 


v w « • * i • * 

+ t-x sin 0^ + y cos e ] 0 - r, 0 (12) 

3 r 6 & 6 * 6 


&U dU 
3r 3 3z 


(13) 


The bracketted quantities ibove correspond to the coordinate OBSDOT 

transformations coded in subroutines XEFIX and YEF1X. These XEF1X 

transforms are used on the true of date satellite velocity YEFIX 
• • 

components x and y. The interested reader should refer to 
Section 3.4 for further information on transformations 
between Earth- fixed and true of date coordinates. 
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I , 
4 
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It should be noted that all quantities in this 
formula, with the exception of those quantities bracket- 
ted, are Earth- fixed values. (The magnitude r is in- 
variant with respect to the coordinate system transforma- 
tions . ) 


The remaining time derivatives are tabulated 

here: 


nocnnT 
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1 

* 
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Right as cension: 


. r. 2 -u 2 r x 

a ■ = — r — 

P U-U3) 


(14J 


Declination: 
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Direction Cosines: t 


p • E-lp 


P 


(16) 


p * N-mp 


m 


P 


(17) 
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c a ~sLLTTE- SATELLITE TRACK*' " 

'•'he f nr, : ..i'-j ta! ■ .itellitc measurement used by 

the GEODYN progiV<-p ’ s sbowr . » Igurc 1. A signal is transmitted 
from a ground tracking sta; -u to one satellite where it is 
then relayed to a second • voilite . The second satellite in 
turn relays the signal back to the fiTst satellite where it is 
relayed to the original ground station. The fundamental measure- 
ment mad'' is the transit time for this relay process. Properly 
corrected for various time delays, this measurement can be trans- 
formed into the sum of the range from the ground station to the 
first satellite and the range from the first satellite to the 
second satellite. The time rate of change of this- measurement 
is also handled by the GEODYN program. 


6.4.1 Satellite-Satellite Tracking Measurement Calculations TKOS 

Given the ephemerides of the two satellites, the range 
sum type measurement can be calculated in a rather straight- 
forward manner. The most important aspect of the calculation 
is to insure that the correct times are used for the satellites 
and ground station. That is, transit times and transponder 
delays must be correctly accounted for. 

To see the times needed for the range sum calculation, 
refer to Figure 1. Let 

R s (t) ■ the range sum measurement at time t 

Rj u ■ the up -link range from the ground to the 

relay satellite 

R£j ■ the relay satellite-tracked satellite range 


I 


R «'«V S»,e!li|» 


( 
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Rp u = the tracked satellite-relay satellite range 

R^ a the down- link lauge from the relay satellite 

to the ground 

R (t) ,R, (t) , (t)= the range vector from the center of the 

earth to the ground station, relay satellite, 
and tracked satellite, respectively, at 
time t 


At Zd 

At 2u 


At 


Id 


* the transponder delay in the relay satellite 

» the transponder delay in the tracked satellit 

* the transit time for the range R lu 

* the transit time for the range 

* the transit time for the range R 2 U 

= the transit time for the range R- . 

w id 


The range sum measurement is expressed in terms of the range 
components as 

" R lu * R 2d * R 2u * R ld (1) 

Each of the ranges on the right hand side is a function of two 
different times. Expressing the ranges in terms of the range 
vectors from the center of the earth and explicitly indicating 
the times, the measurement R_ is expressible as 


/ 


( 

V 


* 


4 

t 



« s (t) - iRjCt-Atjj) - K g (t)| 

+ l lf 2 (t - At ld‘ d r At 2u ) • 

( 2 ) 

+ |R 1 Ct-At ld -u 1 -At 2 u -d 2 -At 2d ) - (t-At ld -d 1 -At 2 u -d 2 ) ! 

♦ |lT 1 Ct*At ld “2d 1 -At^ ,-d 2 -At 2d ) 

-'R g (t-At ld -2d 1 -At 2u -d 2 -At 2d -At lu )| 


This expression shows that the ground station and satellite 
positions must each be known for several different times. 
Summarizing: 


a. Ground station position needed at times 

1. t 

2. t At id" 2d l‘ At 2u‘ <i 2' At 2d’ At lu 

b. Relay satellite position needed at times 


1* * At ld 

2. t - At ld -dj 

3. t - Atjj-dj- At 2u -d 2 -At 2d 

4. t At id‘ 2d i" At 2u’ d 2’ At 2d 


c. Tracked satellite position needed at times 


2 . 


t - At 
t - At 


ld“ d r At 2u 

ld‘ d r At 2u 


d 


2 
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The transponder delay which is most critical is that of the 
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the range rate between the relay and tracked satellite is 
expected to be much higher than the ground-relay satellite 
range rate. This maximum rate can be only on the order of 
5 x 10 3 m/sec, however, and a 4 ysec transponder delay would 
be necessary to introduce a measurement computation error 
of 1 cm. Since actual S-band transponder delays are generally 
no longer than this, we may neglect transponder delays in 
the measurement calculation and still retain accuracies at the 
centimeter level. 


With the neglect of transponder del'ys, we are left 
with 2 times for which the ground station position must be 
computed, 2 times for which the relay satellite position must 
be computed, and 1 time for which the tracked satellite 
position must be computed. Eqn. (2) can then be written in 
the slightly simpler looking form: on> 0C IBIUTY OF THE 

Xnal PAGE IS POOR 


2R 


s (t) - I^Ct-At^)- K g (t)J 
♦ l^2^ t ‘ At ld" At 2u^ ’^1 


TKOSTA 


(3) 


♦ |^ 1 (t-At ld -At 2 u -At 2 d )-T^ (t-At ld -At 2 u -At 2 d *At lu ) | 


This is the form used by GEODYN to calculate the range sum 
measurement. The range sura rate measurement is calculated from 
the time derivative of this expression. To see how this cal- 
culation is performed, note that, e.g. , the final down leg 
range is 

lUjft-Atjp-n col- uRjtt-Atjp-^djj.mjft-Atjp-RjConWi 


and that its time derivative is 


d 

St 


( 4 ) 


J {[RjCt-AtjjJ-R (OI-tRjCtrMjjJ-R (t )]} 1/2 


The calculation thus requires the satellite velocities, and 
the station inertial velocity, at the same times as were 
needed for the range sum computation. The satellite velocities 
are always computed by the GEODYN integrator along with the 
satellite positions, so only the station inertial velocities 
are needed as additional input to the range sum rate calculation. 


6.4.2 Partial Derivative Calculations for Satellite-Satellite 
Tracking Measurements 

Differential corrections for epoch element and force 
model parameter errors require. the computation of the partial 
derivatives of the measurements with respect to these adjusted 
parameters. Let y be one of these parameters. Then, since 
the range and range rate measurements are explicit functions 
of the satellite coordirates only, the partial derivatives 
of R., e.g., can be written from Eqn. (1) as 


aR s e 1 
W 7 


v * 


where 


aR lu 3R 2d aR 2u aR ld 

WJT 1^7 

aR 2d 8R 2ul a *21 

Wl *T“ 


w 
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x 2i * T6 tIie cartesian position coordinates 

of the relay and tracked satellite, respectively. 
Summation over i from 1 to 3 IS implied. 



Eqn. (5) is shown in a somewhat simplified form, since the 
different range sum components depend upon the satellite 
coordinates at slightly different times. For partial derivative 
computations, however, this slight vai..e difference is negligible 
The partial derivatives of the satellite coordinates with 
respect to the y parameters are obtained by independently inte- 
grating the appropriate variational equations for each satellite 
in the same manner in which GEODYM integrated these equations 
for one satellite. 


Eqn . 

C5) 

can be simplified somewhat by noting 

that 


,R 2 d 


m 

3R 2d 

(6a) 


8X li 


8 hi 


8R 2u 
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Using (6a) - 

(6d) , Eqn. (5) can be written 
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TT 

3R ld 3X li . * R 2d ( 9X li - 3X 2i | 

uqi ” + \w~ vr) 
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(7) 


and is a sufficiently accurate for* for the range sun partial 
derivative calculation. 


The partial derivatives of the rirtgc sum trite measure- 
ments are calculated in a similar manner, except that velocity 
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partials must now be included. Thus, if down leg rate partials 
arc approximately equal to up leg rate partials, 
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PCE MEASUREMENTS TYPES 


PREDCT 


The PCE measurement types are sets of elements pre- 
cisely determined in previous GEODYN orbit determination 
runs . 

The inertial Cartesian elements obtained from inter- 
polation of the integrator output are used as the calculated 

• • t 

measurements for PCE types, x,y,z,x,y,z. 

The partials of these measurements are 
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The osculating elements obtained by conversion of the 
above mentioned Cartesian elements are used as the calculated 
measurements for PCE types, a,e,i,0,m,M. 

The partials for these measurements are given in 

Section 11.4. 
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6.6 VLB I MEASUREMENT TYPES 

TNOSTA 

The geometry for the VLBI measurements used by the 
GEODYN program is shown in Figure 1. A signal is transmitted 
from one satellite to two ground stations, 

VLBI Time Delay Measurement Calculation: 



t - is the time delay measurement. 

g 

tj - is the light time for the first ground station. 

* is the light time for the second ground station. 
Pj - is the first station-satellite Tange. 

P 2 - is the second station-satellite range. 


c - is the velocity of light 



Partial Derivative 



X 

-3p 2 ^ 

3p il 

(2) 

8r i 

c 

ar . 

L X 

aTTj 


3 p 2 

and the partials 

3 p i 

, are given in Section 6.2. 
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VLBI Fringe Rate Measurement Calculation: 


V t‘ t ~ ["2 ' Pi] 


where 


f - is transmitter frequency. 

P 2 - is the time derivative of P 2 
Pj - is the time derivative of P^. 

Partial Derivative: 


a Al 

Srj . clarj Sr^J 
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ip 2 

where the partials — are given in Section 6. 2 
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6.7 AVERAGE RANGE RATE MEASUREMENT TYPES 

Figure 1 illustrates the geometry of the average range 
rate measurement t?pes. A signal is transmitted from a 
transmitter to a satellite, then a ground station receives 
the signal from the satellite, and, 


TWO ST A 


P T - is the transmitter-satellite range 

p - is the satellite-receiver range 
K 

K - is the ppsition vector of the receiver 

£ . is the position vector of the transmitter 

£ _ i S the position vector of the satellite. 

If t is the recorded time of the end of the doppler 
counting interval at the receiver and, if T is the length 
of the counting interval, then the average range rate 
measurement is 


P .■ 


P R (t 6 ,t 5 ) ♦ P T (t s ,t 4 ) - P R (* 3 ,t 2 ) - P T (t 2* t l ) 


2T 


( 1 ) 


Where it is necessary to iterate for the satellite 
and transmitter times, 


*5 “ *6 


p R (t 6» t 5 ) 
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Figure 1. Geometry for Average Range Rate Measurement 
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p R^ t 3 ,t 2 ' ) 
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p T( t 2 ,t l^ 


and where 


p R^ t 6 » t 5 ^ ~ 




p T ( t 5 » t 4) = l^xC t 4) " "sC^I 

p R^ t 3 ,t 2^ “ ^R^3^ ' ^S^2^ 

P T (t 2 .t 1 ) = - Eg(t 2 )l 


A two-way average range rate measurement is a special 
case of the three-way average range rate measurement (i.e., 
the receiver and the transmitter are the same). Therefore, 


P T • P R , «r 
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The Partial Derivatives are 
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SECTION 7.0 
data PREPROCESSING 


The function of data preprocessing is to convert 
and correct the data. These corrections and conversions 
relate the data to the physical model and to the co- 
ordinate and time reference systems used in GEODYN 
The data corrections and conversions implement'd in 
GEODYN are to 

• transform all observation times to A1 time 
at the satellite 

• refer right ascension and declination ob- 
servations to the true equator and equinox 
of date. 

• correct range measurements for transponder 

delay and gating effects 

• correct SAC right ascension and declination 
observations for diurnal aberration 

• correct for refraction 

• convert TRANET Doppler observations into 

range rate measurements. 

4 

These conversions and corrections are applied to the data 
on the first iteration of each arc. Each of these pre- 
processing items is considered in greater detail in the 
subsections which follow. 
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7.1 


TIME PREPROCESSING 


The time reference system used to specify the 
time of each observation is determined by a time 
identifier on the data record. This identifier also 
specifies whether the time recorded was the time at 
the satellite or at the observ ; ng station. 

The preprocessing in GEODYN transforms all 
observations to A1 time in either GEOSRD or DODSRD. 

If the time recorded is the time at the station, it 
is converted to time at the satellite. This con- 
version is applied in subroutine PROCES using a cor- 
rection equal to the propagation time between the 
spacecraft and the observing station. The station- 
satellite distance used for this correction is computed 
from the initial estimate of the trajectory. 

There is special preprocessing for right 
escension and declination measurements from the GEOS 
satellites when input in National Space Science Data 
Center format. If the observation is passive, the 
image r^orded is an observation of light reflected 
from the satellite and the times are adjusted for 
uropagation delay as above. If the observation is 
active, the image recorded is an observation of light 
transmitted from the optical beacon on the satellite. 
The beacons on the GEOS satellites are prograauned to 
produce a sequence of seven flashes at four second 
intervals starting on an even minute. For the active 
observations, the times are set equal to the programmed 
flash time with a correction applied for known clock 
errors (Reference 1), plus half a billisecond, the time 
allowed for flash buildup. 
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DODSRD 

GEOSRD 

PROCES 


GEOSRD 




f 

The corrections for the active observations are 
applied in GEOSRD, which calls SATCLC and SATCL2 to 
evaluate the corrections for GEOS 1 and GEOS 2, re- 
spectively. These routines compute the correction by 
simple linear interpolation in a table of known errors 
in the satellite on-board clock. 
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REFERENCE SYSTEM CONVERSION TO TRUE OF DATE 


DODSRD 


The camera observations, right ascension and 
declination may be input referred to the mean equctor 
and equinox of date, to the true equator and equinox 
of date, or to the mean equator and equinox of some 
standard epoch. The GEODYN system transforms these 
observations to the true equator and equinox of date 
in subroutines GEOSRD and DODSRD. The necessary 
precession and nutation is performed by subroutine 
EQUATR. 


7.3 TRANSPONDER "ELAY AND GATING EFFECTS 

The range observations may be corrected for 
transponder delay or gating errors. Tf requested, the 
GEODYN subroutine PROCES applies the corrections. 

The transponder delay correction is computed as 
a polynomial in the range rate: 


Ap « a Q ♦ a x p + a 2 (p)“ 


where a^, a^ , and a 2 depend on the characteristics of 
the particular satellite. 


A gating error is due to the fact that actual 
range measurements are either time delays between 
transmitted and receive, radar pulses or the phase 
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shifts in the modulation of a received signal with 
respect to a coherent transmitted signal. Thus there 
1 r» rnrrprt iv identifying the Te- 

is me iu a x a - - * 

turned pulse or the number of integral phase shifts. 

The difference between the observed range and the computed 
range on the first iteration of the arc is used to deter- 
mine the appropriate correction. The correction is such 
that there is less than half a gate, where the gate is the 
range equivalent of the pulse spacing or phase shift. The 
appropriate gate of course depends on the particular 

station. 


7.4 ABERRATION 


PROCES 


Optical measurements may require corrections (Refer- 
ence 2) for the effects of annual aberration and diurnal 

aberration. 


Annual Aberration 

The corrections to right ascension and declination 
measurements for annual aberration effects are given by 



20V 5 (cos a’ cos * cos ^ + sin a* sin *) 


cos 6’ 



4 • 6* - 20VS (cos • ccs e T (tan c T cos 6' -sin a' sin 6 r ) 
♦ cos o' sin 6' sin •] 
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where 
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o - true right ascension of the satellite 


a* - observer right ascension of the satellite 


6 - true declination of the satellite 


6* - observed declination of the satellite 


e.j. - true obliquity of date 


• - geocentric longitude of the sun in the ecliptic 

plane 


Diurnal Aberration 


C 


The corrections to right ascension and declination 
measurements for diurnal aberration effects are given by 


a ■ o' + 07320 cos cos h. sec 6* 

9 


5 - 5' ♦ 07320 cos sin h sin 6’ 
• 21 


where 


4> ' - geocentric latitude of the statibn 


h - local hour angle measured in the westward 
s 

direction from the station to the satellite 


a - true right ascension of the satellite 


o' - observed right ascension of the satellite 
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6 - true declination of the satellite 

6’ - observed declination of the satellite 


7.5 REFRACTION CORRECTIONS 

The GEODYN system can apply corrections to all 
of the observational types significantly affected by 
refraction. The corrections requested are applied by 
subroutine PROCES. 


Right Ascension and Declination: 

Optical measurements may require corrections 
(References 3, 4, 5) for the effects of parallactic re- 
fraction. These corrections are given by 


a * o’ - AR sin q/cos l 
6*6* - AR cos q 

where the change in the zenith angle, AR> in radians is 
given by 


AR 


0.435 (4.84813) tan Z Q 
p cos Z Q 


(1 . e (-l.S8S) ID'* 


and 


a - true right ascension of the satellite 
a' - observed right ascension of the satellite 
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cos Zoj 
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true declination of the satellite 


PROCES 
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6’ - obseived declination of the satellite 


Z Q - observed zenith angle in radians 

p - range from the station : - the satellite in 
meters 


i 

i 


i 

i 

! 
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q - parallactic angle in radians 

The parallactic angle q is defined by the intersection 
of two planes represented ty their normal vectors F^ and F 2 . 


A A 

* C p x u 

A 

f 2 * v x 11 


where 


A 



(0,0,1) 


v - unit local vertical at the station 



u - unit vector pointing from the station to the 
satellite in inertial space. 


t* 

(4 
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Therefore, the sine and cosine of the parallactic angle 
are given by 

* 

A A 

cos q * Pj . ?2 

A A 

sin q * Pj . P 2 
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Pi - unit vector in the I 7 ! direction 
P 2 * unit vector in the ^ direction 

and 

A 

~ P, x u 

I p = _± 

; 3 M’j x u | 

I 

t 

The parallactic angle, q, is measured in the clockwise 
direction about the station-satellite vector (i.e.*, a left- 
handed system is used to define this angle). All vectors 
and vector cross products used in this formulation conform 
to a right-handed system. 


( 


Range : 


PROCES 
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The refraction correction applied to CNES laser 
range data is 


A p n 

A p * . - — ““T 

sin E t + (cot E t ) 10 

and the correction applied to range data from all other 
tracking systems is 



2.77n 

Ap - 

328.5(0.026+sin E^) 


( 4 ) 
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where 

A p - is that correction associated with a range 
n 

observation measured along the direction of 
the satellite zenith, and is provided along 
with each observation on the data tape, 

E^ is the elevation angle computed from the 
initial estimate of the trajectory 


PPM deviation fx„m unity of the surface 
index of refraction; if this value is not 
specified, it is assumed to be 328.5. 


Range Rate: 

■ For range-rate , the correction Ap is derived from 
the range correction: 


Ap 


2.77n,. cos 
. 328 .S(0.026+sin E^) 


7 


(5) 
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where 


is the computed rate of change of elevation. 


Elevation: 

For elevation observations the correction AE^ 
as computed as follows: 

n 10* 

AE, « 2 (6) 

* 16.44+930 tan E^ 

Azimuth is not affected by refraction. 

Direction Cosines: 

The corrections At and Am are derived from the 
elevation correction: 


At « 

-sin A z sin 

(E t ) 


(7) 

Am ■ 

-cos A z sin 

(E t ) 


(8) 
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where A z is the azimuth angle computed from the initial 
estimate of the trajectory. 


PROCES 


X and Y Angles: 

For X and Y angles the corrections AX and AY are 
computed as follows: 


sin A-AE„ 

AX * Z —± 

(sin 2 + sin 2 A^cos 2 E^) 



cos A z sin E £ AE^ 
y i-cos 2 A 2 cos 2 E^ 


( 9 ) 


( 10 ) 


Note that these are . Iso derived from the elevation correction. 


7.6 7RANET DOPPLER OBSERVATIONS 

TRANET Doppler observations are received as a GEOSRD 

series of measured frequencies with an associated base 
frequency for each station pass. Using the following 
relationship, the GEODYN system converts these observa- 
tions to range rate meas- rements in subroutine GEOSRD: 


P 


c < F lf F M> 


M 


U) 
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Fjj is the measured frequency, 

F R is the base frequency. 


and 


c 


is the velocity of light. 


SATELL.TE -SATELLITE TRACKING DATA PREPROCESSING 


TffOM A 
UP UCh. 


The preprocessing on the satellite-sat >11 ite tracking 
involves the determination of all the appropriate transit 
times. Because of the station -satellite and inter -satellite 
dis tances, this process must be performed iteratively. The 
required times are computed during the first iteration and 
are then stored for use in subsequent iterations. 

The satellite-satellite tracking measurements are a*so 
corrected for tropospheric refraction. The corrections made 
here are identical to those which would be made on range and 
range rate measurements to the relay satellite only. Although 
it is theoretically possible for signals from the relay to 
low altitude satellite to pass through the atmosphere, such 
tracking would occur at reduced signal intensity and v.'Ould 
be equivalent to the low elevation tracking of satellite from 
ground based stations. Such date is seldom used in orbit 
estimation. 

The standard procedure for transponder delay 
corrections on satellite-satellite tracking is to use 
block data constants for each satellite, with a satellite 
ID used to identify the appropriate block data entiie*. 

Since constants for the transponders to be used for 
satellite-satellite tracking are not presently available 
the block data entries must be modified appropriately v/hen 
the data becomes available. 
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SECTION 8.U 

FORCE MODEL AND VARIATIONAL EQUATIONS 


A fundamental part of the GEODYN system requires 
computing positions and velocities of the spacecraft 
at each observation time. The dynamics of the situa- 
tion are expressed by the equations of motio. - , which 
provide a relationship between the orbital elements 
at any given instant and the initial conditions of' 
epoch. There is an additional requirement for varia- 
tional partials, which are the partial derivatives of 
the instantaneous orbital elements with respect to the 
parametc.s at epoch. These partials are generated 
using the variational equations, which are analogous 
to the equations of motion. 


8.1 EQUATIONS OF MOTION 

In a geocentric inertial rectangular coordinate 
system, the equations of motior for a spacecraft are of 
the form . 

.. ur 

r - - -r ♦ * (1) 

r 

where 

r is the position vector of the 
satellite . 
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y is GM, where G is the gravitational constant 
and M is the mass of the Earth. 

"K is the acceleration caused by the 
asphericity of the Earth, extra- 
terrestrial gravitational farces, atmos- 
pheric drag, and solar radiation. 

This provides a system of second order differential 
equations which, given the epoch position and velocity com- 
ponents, may be integrated to obtain the position and velocity 
at any other time. This direct integration of those 
accelerations in Cartesian coordinates is known as. 

Cowell's method and is the technique used in GEODYN's 
orbit generator. This method was selected for its 
simplicity and its capacity for easily incorporating 
additional perturbative forces. 

There is an alternative way of expressing the F 

above equations of motion: 


r * V U * A p ♦ ^ 


S. 


where 


( 2 ) 




U is the potential field due to gravity, 

Xp contains the accelerations due to drag, 
and 

contains the accelerations due to solar 
radiation pressure. 

8 . 1-2 
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1 


r 


I 


This is, of course, just a regrouping of terms coupled 
with a recognition of the existence of a potential field. 
This is the form used in GEODYN. 

The inertial coordinate system in which these 
equations of motion are integrated in GEODYN is that 
system corresponding to the true of date system of ol*o 
of the reference day. The complete definitions fcr these 
coordinate systems (and the Earth- fixed system) are 
presented in Section 3.0. 

•« 

The evaluation of the accelerations for r is 
controlled by subroutine F. This evaluation is performed 
in the true of date system. Thus there is a requirement 
that the inertial position and velocity output from the 
integrator be transformed to the true of date system for 
the evaluation of the accelerations, and a requirement to 
transform the computed accelerations from the true of date 
system tq the inertial system. These transformations are 
performed by subroutine REFCOR (which controls the pre- 
cession and nutation routines, PRECES and NUTATF.) and is 
controlled by subroutine F. 


8.2 THE VARIATIONAL EQUATIONS 

The variational equations have the same relationship 
to the variational partials as the s-.cellite position vector 
does to the equations of motion. The variational partials 
are defined as the t where x t spans the true of date 

position and velocity of the satellite at a given time; i.e. t 


F 

REFCOR 


VEVAL 
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f x . y , z . x , y , z } ; 


VEVAL 


Y 

' t 


and g spans the epoch parameters; i .e., 


W z o 


the satellite position vector at 
epoch 


W z o 


the satellite velocity vector at 
epoch 


the satellite drag factor 

the time rate of change of the drag 
factor 

the satellite emissivity factor 


C ,S 
nm' nm 


gravitational harmonic coefficients 
for each n, m pair being determined 


X. 


surface density coefficients 


Let us first realize that the variational partials 
may be partitioned according to the satellite position 
and velocity vectors at the given time. Thus the re- 
quired partials are 


dr 3r 
d? * d? 


( 1 ) 
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r is the satellite position vector (x,y,z) 
in the true of date system, and 

• • • * 
r is the satellite velocity vector (x,y,z) 

in the same system. 


The first of these, — , can be obtained by the double 
integration of 


© 


or rather, since the order of differentiation may be 
exchanged, 


Note that the second set of partials, , may be obtained 

SB 

oy a first order integration of — . Hence we recognize 

^ Sr 

that the quantity to be integrated is — . Using the second 

form given for the equations of motion in the previous 
subsection, the variational equations are given by 


■M 0 
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3r 3 


( 4 ) 


VEVAL 


30 30 


(VU 


'R 


+ V 


where 

U is the potential field due to gravitational 
effects 

is the acceleration due to radiation pressure 
is the acceleration due to drag 

The similarity to the equations of motion is now obvious. 


At this point we roust consider a few items: 


T 

* V .Vu 


• The potential field is a function only of 
position. Thus we have 



• The partials of solar radiation pressure 
with respect to the geopotential co- 
efficients, the drag coefficient, and the 
satellite velocity are zero, and the par- 
tials, with respect to satellite position, 
are negligible. 

• Drag is a function of position, velocity, 
and the drag coefficients. The partials, 
with respect to the geopotential coefficients 
and satellite emissivity, are zerc, but we 
have 

a*D ^ 3A D 3x t + * C D + 

3ir 3x t 31 3 Cjj 3? 3C D 3F 


( 6 ) 


f 

Let us write our variational equations in matrix 
notation. We define 


c 




c 

I 


n to be the number of epoch parameters in 8 

F is a 3 x n matrix whose j** 1 column vectors 

are 3r 

FFT 


U 


2c 


is a 3 x 6 matrix whose last 3 columns are 
zero and whose first 3 columns are such 
that the i, element is given by 

a 2 u _ 

57 T 5f j 


m 


is a 3 x 6 matrix whose i ^ column is defined 

. **n 

by D 

dx+ 

is a 6 x n_matrix whose x tow is 

given by . Note that X_ contains the 
TB7 

variational partials. 


m 


f ' is a 3 x n matrix who. e first six columns 
are zero and whose last n-6 columns are 
such that the i, j th element i‘ given by 

. ■y|~ (VU ♦ . Note that the first six 

columns correspond to the first six elements 
of F which are the epoch position and velocity. 
(This matrix contains the direct parV/tls of 
x t with respect to F. ) 


8 . 2-6 


f 


I 


i 


n 


VEVAL 


1 ( 

We may now write 


i 


REPRODUCIBILITY OF THE 
ORIGINAL PAGE IS POOR 


i , F ■ l“2c * D r> ™ 

! 

] 

! i 

This is a matrix form of the variational equations. 

i 

i 

i 

i Note that U 2c , D r , and f are evaluated at the 

j current time, whereas X m is the output of the integra- 

i tion. Initially, the first six columns of X m plus 

! the six vows form an identity matrix; the rest of the 

; matrix is zeroffor i»j,X m *1; for i?* j , X m *0) . 

Because each force enters into he variational 
equations in a manner which depends directly on its 
.. - model, the specific contribution of each force is dis- 

cussed in the section with the force model. We shall, 
however, note a few clerical details here. 



The task cf computing these variational equations 
in the GEODYN system is largely accomplished by sub- 
routine VEVAL. The matrix dimensions given are for 
notational convenience; empty rows and columns are not 
programmed . 


The above equation is also applied in subroutine 

\ PRELCT to "chain the partials back to epoch," that is, to 

relate the partials at the time of each set of measure- 

| ments back to epoch. 

• \ * 
sr 



VEVAL 


PREDCT 


August 33 , 3 973 


3 x 

The matrix for , X m above , is initialled in C-RiilT 


subroutine ORBIT. 


D71 


.... . _ r _..t nvi DRAG D650 

Tne contriDULiuns uj. auux 

EGRAV , F, SURDFN, and RF.SPAR will be discussed as part PRAG 

of the following subsections. The matrices U 2c and f will f 
be referred to in each subsection as though the particular RESPAR 
force being discussed had the only con.rirmtion, SURDF.N 


( 


0 
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THE EARTH'S POTENTIAL 


6.7 

In GEODfN che Earth's potential is described by 
the combination of a spherical harmonic expansion and a 
surface density layer. Generally, however, the spherical 
harmonic expansion is used exclusively and no surface 
density terms are included. 



8.3,1 Sph erical Harmonic Expansion 

The Earth's potential is most conveniently ex- EGRAV 

pressed in a spherical coordinate system as is shown 
in Figure 1. by inspection: 

• <t>', the geocentric latitude, is the angle 
measured from OX, the projection of OF in 
the X-Y plane, to the vector OP. 

• A, the east longitude, is the angle measured 
from the positive direction of the X axis 

to OX. 

• r is the magnitude of the vector OF. 

Let us consider the point P to be the satellite E'lRAV 

position. Thus, OF is the geocentric Earth-fixed satellite 
vector corresponding to r, the true of date satellite 
vector, whose components are (x,y, 2 ). The relationship 
between the spherical coordinates (Earth- fixed) and the 
satellite position coordinates (true of date) is then 
given by 



£.3-2 

i 





EGRAV 


where 6 is the rotation angle between the true of date 

g 

system and the E^rch-fixed system (see Section 3.4). 


The Earth's gravity field is represented by the 
normal potential of an ellipsoid of revolution and 
small irregular variations, expressed by a sum of 
spherical harmonics. This formulation, used in the 
GEODi'N system, is 


nmax n 


U = f l ♦ 
r 


e m * h >) [‘ 


n=2 m=-0 


cos mA + S sin mX 

nin nm 


( 4 ) 


where 


G is the universal gravitational ccnstar.t. 


M is the mass of the earth, 

r is the geocentric satellite distance, 

nmax is the upper limit for the summation (highest degree), 


a e is the Earth's mean equatorial radius, 
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n 



is the satellite geocentric latitude, 


EGRAV 


X is the satellite east longitude, 


P™(sin indicate the associated Legendre 


functions, and 


(! and are the 
nr. nm 


gravitational 


coefficients . 


The relationship., between th^ normalized co- 
efficients „) *nd the denormalized coefficients 

are as follows: 


DENORM 


[rn-m)!(2n*lK2-« om ) j 1/2 

(n+m) ! 


where 


is the Kronecker delta, 
om ' 


6 = 1 for m=0 and 6 =0 for myO. 

om om 


A similar expression is valid for the relationship 

between 3? and S_ m . This conversion factor is com- 
run nm 

puted by the GEODYN system function DENORM. 
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\ 


( 


( 
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The gravitational acceleration* in true of date co- 
ordinates (x,y,z) are computed f . >m the geopotential, 
U(r,<t>\ X) , by the chain rule; e.g., 


3 U 

dr 3U 

36 

3U 

3X 

— 

+ 

+ 

— 

3r 

dX 3$' 

3x 

3 X 

3x 


( 6 ) 


The accelerations y and z are determined likewise. The 
partial derivatives of U with respect to r, and X are 
given by 


3U 

3r 


GM 

7 


1 + 


nmax , n 
a e\ 



n=2 


\r, 


]C (C nm 


cos mX 


(7) 


m -*0 


+ S nm sin mX) (n + 1) P“ (sin <p*^ 


3U GM 
3X 


nw nmax . n n 

“E 0 E 

n*2 ' r/ m*0 


(S cos raX - C, m sin mX) (8) 

v nm * v 3 


nm 


m P® (sin 40 


nmax , . n n 


3U GM /a\** " 

7, 1 - L H £ < C n* cos * X + S iim sin BX) ' 5) 

»♦' r £i Xrl 


P n +1 (sir * *’ 5 * m tart ♦' ^ (sln ft 
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The partial derivatives of r, and X with respect to 
the true of date satellite position components are 


3r 


r i 


3 


r 


1 



ax 

9r. 


- 

1 

‘ r i ♦ 

32 j 

' ~T~ * 

a ! 

r 

3r i J 

r Sy 

y 3x 1 

t ar. 

x ar. 1 


( 10 ) 


(ID 


EGRAV 


,n — T 


( 12 ) 


The Legendre functions are computed via recursion 
formulae : 

Zonal s: m*0 

P° (sin 40 = - [(2n-i) sin A ' 1 

n n L 

(n-1) P®. 2 (sin ♦oj 
P° (sin 3 sin 4>' 

Tesserals: m^O and m<n 

P m (sin *0 « p „ ? ( sin ^ + (2n ‘ lj cos *' P n-1 (s 
n n- l 

pj (sin 40 - COS 
Sectorals: m**n 

P m - (2n-l) cos ♦' Pjj’.t (sin 40 



The derivative relationship is given by 


EGRAV 


d /p® (sin fn - P^[ +1 (sin <*0 - m tan ♦' Pj/.sin 40 

d*' ' n 11 (18) 


It should also be noted that multiple angle 
formulas are used for evaluating the sine and cosine 
of mX. 


ECRAV 

VEVAL 


These accelerations on the spacecraft are com- 
puted in subroutine EGRAV. Arrays containing certain 
intermediate data are passed to subroutine VEVAL for 
use in the computations for the variational equations. 
These contain the values foT: 


GM 

r 



( 19 ) 


P* (sin ♦*) 
sin mX 
cos mX 
m tan *' 


for each » and/or n. 


VEVAL 


The following discussion relates primarily to 
the mathematical formulations utilized in subroutine 
VEVAL . 


The variational equations .*equiie the computation 
of the matrix U > whose elements are given by 


( 


U 


2c 



2 

3* U 


3r. ar. 
i J 


( 20 ) 


where 

u (x, y, z), the true of date satellite position. 
U is the gc opotential . 


Because the Earth's field is in terms of r, sin 
and A, we write 

3 

U 2c “ C 1 U 2 C 1 + 

k-1 


3U 


3e n 


C 2k 


( 21 ) 


where 

e^ range-' over the elements r, sin and X 

U, is the matrix whose i, j tn element is given 

by ** U 
^e^TeJ 
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C, is tne mull ix WuCi 

1 . 3e- 

by i 


is aiven 


is a set of three matrices whose i, j 


elements are given by 3 

5rT SrT 

*■ J 


We compute the second partial derivatives of the 
potential U with respect to r, <J>\ and X: 


3 2 U 2GM GM 


n n 


GM GM /» e Y " 

7*7 S (n+1) ,n * 2) ( 7 ) £ 


(C cos + ^nm s * n ^ 


run ax m \* 1 »* 


n n 


3r3$' 


(C _ cos mX 
run 


♦ S„ M sin mX) — (p” (sin tf)) 


GM 

* ^ 

nmax 


ya \ n n 

L 

(n+1) 

ft E 


n-2 


rr a«*o 

(-C 

1 nm 

sin mX 

♦ S 

nm 

cos mX) 


8 . S-ll 


1 


3 2 U 

17 


ru nmax Ja .n n 

■"se £« 


n=2 


VPVAT. 


cos mX « S__ sin mX) • 
nm nm 


m*o 


(25) 


3 

3* 


(p® (sin <|f)J 


3 2 U 

d^dX 


noax «, > n n 


CM « /a A* — 

7 r (7) e - ( - c » sin bx 


n-2 


m*o 


(26) 


+ 3 nm cos 


s mX) — (p® (sin ♦')) 
3*’ ' c ' 


3 2 U 

II 7 


nmax » ja n 
n i, 2 ' m*o 


^ nmax , .n n 
GM (&A 7 

m (C cos mX 

nm 


(27) 


♦ S__ sin mX) P® (sin 40 
nm n 


where 



P®* 1 (sin ♦’) - m tan ♦' (sin ♦*) 

" (28) 
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f 


3 ** \ 

~2 ^P® (sin sm - P** +2 (sin - (m+1) tan <*>' P® + “(sin $*) 

- m tan 4> ? |pj| + ^ (sin $0 - m tan $' P^(sin (frOj 

- m sec 2 4> P* (sin ») (29) 


The elements of U 2 have almost been computed. 
What remains is to transform from (r, X) to 
(r, sin This affects «**ly the partials involving 

4 >*: 


( 


3U 

3U 

e«b' 





3 s in <J>’ 

3<t»' 

3 sin 





3 2 U 


H' I 

f3 2 u\ 

3$’ 

3U 

3 2 $' 


3 sin 4>‘ 


3 sin 4? 


(30) 


3 sin <j>' 3<J>' C sin 4>' 


TZ 


(31) 


where 


i 


?♦* 

3 sin <j>' 


sec 4> l 


(32) 



3 ♦' 


3 sin $ 


7 ’ 


sin ♦’ sec' 


(33) 


VEVAL 
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For the and C2^ matrices, the partials of r, 
sin 4>', and X are obtained from the usual formulas: 


r 


f~2 — 2 — 7 

\V+ y 4 +z* 


( 34 ) 



( 35 ) 


( 36 ) 


We have for C^: 



r. 

l 

r 


3 sin <J>* 



- z r . 
i 



1 

T 



3X 1 

= ~1 — 1 

3r. x +y 4 




( 37 ) 


( 38 ) 


( 38 ) 
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The C ov are symmetric. The necessary elements 
are given by 


? 2 r 


3 r. 3 r. 
i J 


ILL i . I !fi 

r r 3r^ 


(39) 


3 sin 

3r. 3r. 
i 3 


3 2 X 


3r. 3 r. 
i J 


3z r f r. 

— i~ 


-2r. 

— 2 — / 
(x^YT 


3 z 


3 z 


3 r. 


r. — 
J 3 


+ r + Z : 

1 3 J 


3y 


X - 


3r i 


3x 


3r. 


(40) 

(41) 


1 

I 

x +y 


3x 3y 3y 3x 


,r j s h 


3 3 r ^ 


If gravitational constants, C n# or are being 
estimated, we require their partials in the f matrix 
for the variational equations computations. These 
partials are 



3 


3C 


nm 




■rO 


sin 


cos (mX) P®(sin 4?) (42) 


(mX) P® (sin ♦*) (43) 


VEVAL 


i 

J 


RESPAR 
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3 


r 


one 

1\JU 


( 


( 


r 


i 



GM /aA" 


COS 


(mX) 


p m+l 

n 


(sin <b') 


- m tan (sin f) | 

i 


(44) 


The partials for S nm are identical with cos (mX) re- 
placed by sin (mX) and with sin (mX) replaced by 
-cos (mX) . 

These partials are converted to inertial true of 
date coordinates using the chain rule; e.g.. 



This particular set of computations is performed by 
subroutine RESPAR. The items which EGRAV computes for 
VEVAL are also available to RSSP.\R and are therefore 
utilized. 
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8.3.2 


Surface Density Layers 


The representation of the earth's gravitational field 
in terms of a simple density layer spread over the surface 
of the earth was first introduced by Koch [Reference 10] in 
1968. Attempts at determining numerical values for surface 
densities on a global scale have been made using both optical 
[Reference 6] and Doppler [Reference 7] data. In some cases, 
the surface densities have been estimated as alternatives to 
the spherical harmonic expansion, and in other cases the 
surface densities are a supplementary contribution to a set 
of "known*' low degree and order spherical harmonic coefficients. 

The surface densities implemented in the GEODYN program 
are basically in the nature of a supplementary potential con- 
tribution. The spherical harmonic field is retained for repre- 
senting the geopotential on a global scale and the surface 
densities can be introduced on either a local or global scale 
into any number of olocks of constant density. That is, the 
fineness of representation of the geopotential via surface 
densities is arbitrarily small, consistent with computer core 
availability and the existence of data for actually resolving 
a large number of surface densities. In addition, t..e capa- 
bility now exists in the GEODYN program for simultaneously 
adjusting both spherical harmonic coefficients and surface 
layer densities. No investigator has apparently yet attempted 
this. When actually making simultaneous adjustments, the 
results must be very carefully interpreted. This problem is 
considered further below in the discussion of constraints. 


8. 3. 2.1 Mathematical Representation of Surface Densities. 
The total potential of the earth 'W* can b-s, somewhat 
arbitrarily, divided into a spherical harmonic*. ^art 'U* and 
a remainder ' T* to be expressed in some other form 


W » U + T 


d: 


with 


(sin <j>) 

+ 1/2 ij’ r ^ cos^ <J> (2) 

where r is the distance from the point of interest to the center 
of mass of the earth and $ ar.d X are geocentric latitude and 
longitude. The last term in ( 2 ) is omitted if the potential 

t 

is beine computed outsid' the surface of the earth. In GEODYN. 
the ma' r gree spherical harmonic coefficient is basically 
arbitr<-r; „ /.urmally being limited to the maximum degree for 
which coefficients are available. 

The potential T can be represented as that of a 
simple layer distributed over the surface of -he earth. Mathe- 
matically, T is then given by the surface integral 

T * / / X d£/ * (3) 

S 

where *• is the distance from a point on the surface to the point 
at which the potential is to be computed, dE is the element 
of surface area, x is the surface density (in units of kg/ m 2 
multiplied by G) , and S is the surface of the earth. Figure 1 
shows the geometry and a portion of the surface areas. To 
numerically evaluate the integral in (3), it is necessary to 
divide the entire surface into blocks of constant density. 

If there are M such blocks, then (S' can be written 
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A 



CC 

JJ 

AE . 


AV f C 


f d) 


where x< is new the average density on the i ' th block and 
the integral is to be taten over the area of the- i ' tir block. 


The integral in Eqn. (4) must be evaluated numerica 1 ly . 
It is evaluated in GEODYN by dividing the area up into four 
blocks of equal area and taking the kernal, 1/8, to be constant 
over earn of these sub-blocks. This is the division which 
has been most commonly used for surface density layers and 
has been shown by Koch (Reference 8 ] to be a quite good 
approximation, genet ally accurate to vnthir. a few percent. 

Rest. . of numerical tests are also given below. 


With the division into sub-blocks, the potential due to 
surface densities is 


T 


M 


E 

i= 1 


4 


E 

3*1 


AE. ./£• - 
ij 


(5) 


where AE^ is the area of the j * th sub-division of the i * th 
block and 8^ is the distance from the center of this sub- 
division to the point where the potential is to be evaluated. 
The acceleration produced by the surface density potential is 
obtained by taking its gradient, 


a) 


surface densities 


M 

*r= £ 

i= 1 


4 

E 

3-1 


4E ij 


( 6 ) 

SURDEN 


The forcing function for integrating the variation equations 
to obtain the sensitivity of satellite position to a particular 
surface density block is obtaining by differentiating Eqn. (6) 
with respect to X- , 

i 


8.J-2U 



(7) 

SURDE.N 


Note that these forcing functions must be computed as part of 
the computation of the surface density acceleration contribution. 
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8. 3. 2. 2 Surface Height Computation. A number of potential 
choices are available for locating the surfaces on which the 
surface densities are to be spread. Such surfaces include 
the spheroid, the geoid, aid the physical surface of the 
earth. The method which has been implemented in GEODYN is 
to locate the density layers on the geoid defined by the 
earth and geopotential model being used in the program. 

The model presently being employed is the SAO 1969 Standard 
Earth [Reference 9] . 

The geoid choice for locating the surface densities is 
the most natural for use in estimating surface density values 
in blocks restricted to ocean areas, as might be one of the 
initial uses of the GEOS-C altimeter data. For complete global 
density layers, and perhaps incorporating measurements of 
surface gravity, some other surface may be more convenient. 


8. 3. 2. 3 Layer Model Quadrature Errors. The 
piocess of approximating the integral over the area 
of a surface density block by 4 sub-blocks with 
the kernal estimated at the center introduces some 
error into the integration of surface density effects 
on the orbit. Koch [Reference 3] has investigated the 
error introduced by dividing tne blocks into only 4 sub- 
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blocks, and concluded that, errors generally less than a few 
percent were introduced. 

A test was made in GEODYN to determine the effects of 
different divisions of a 20° x 20° block for a satellite 
of 500 nm altitude passing directly over the center the 
block. The results for a subdivision into 4, 9 and 16 blocks 
are shown in Figure 2. This Figure shows that the 4-block 
subdivision does indeed introduce substantial error, but only 
when the satellite is directly over the center of the block. 

It should be noted that a 29° block size is much !irger than 
would normally be considered for the fine detail representation 
of the geopotential. A division into 20° x 20° tlccks on a 
global scale is, of course, a reasonable possibility. 

Figure 3 shows the acceleration effects due to the 20° 
x 20° surface density layer for a complete revolution of the 
500 nm satellite. It will bo noted that the effects are quite 
localized, as is indeed one of the advantages of the surface 
density representation. There is a large perturbation then 
the satellite is directly cveT the block. There is a definite 
but rather small perturbation when the satellite comes over 
the next revolution about 10° away from the edge of the block. 
Otherwise, tne effects of the blocks are rather constant and 
small . 
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S.3.2.4 Constraints. For several reasons, it is necessary 
to apply certain constraints to the surface density adjust- 
ments which are to be allowed. That this is necessary can 
be seen by noting that the total surface density potential 
can be <xpressed in terms of a spherical harmonic series. 
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which is of the identical form as the global spherical 
harmonic expension given by Eqn. (2), except that the 
expansion is now infinite. It is most significant, how- 
ever, that the surface density expansion could actually 
include the equivalent perturbations of the normal 
spherical harmonic set of coefficients, and that both 
numerical and interpretation problems can arise if both 
spherical harmonic coefficients and surface dcnsitites 
are allowed to adjust simultaneously. 

It may also be noted that fir^' degree coefficients 
in (8) would not, in general, be zero. It is thus 
necessary to force the distribution of densitites to be 
such that these coefficients are i . o in order to avoid 
moving the center of mass of the earth. 

The form which constraints should take can be 
found by expressing 1/i in Eqn. (4) in terms of spher- 
ical harmonics ar.d identifying coefficients P® (sin $) 

cos mX and P® (sin 4) sin mX. Schwarz (Reference 15] has 
n 

shown that this leads to expressions for C' „ a:*d S' ^ of 

nm nm 
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where f s ! if m=0 and zero otherwise. This set of 

om 

integrals can be obtain.^ numerically by breaking the 
area up into sub-blocks as was done for the acceler- 
ation calculation. 

The constraint equations are obtained by setting 
C' „ and S’ _ equal to zero for every spherical harmonic 
coefficient to which the surface densities should not 
contribute. In GEODYN, the default set of zero coef- 
ficients has been set to C'^q, s H' Additional 

constraints (as, c,g., no contribution to 8th degree 
or lower degree coefficients) can be imposed upon 
.input option. 

The GEODYN implementation of constraints is through 
the solution for a number of densities equal to the total 
number of densities adjusted less the number of constraint, 
equations. The normal matrix thus contains only inde- 
pendent densities and core requirements are minimized. 

The procedure for eliminating densitier* is seen by 
writing the constraint equations obtained from (y) as 
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for m n 
n < N' 


where N' is the maximum degree coefficient unaffected by 
the surface density layers. 

The set of Eqns. (10) can be written formally 
as 



j = 


1 , 


M r 


where the are given by the surface integrals in (10) 
and M* is the number of constrain': equations. The 
number M' is related to N' by 


M' * N'(N'+2), 


(10a) 


(10b) 
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as follows from the number of C' _ and S’ _ coefficients 

nm nm 

for which n £ N' and which are not identically zero. On 
the assumption that ri f M, (ii) can be written 



i-1 i-M‘+l 


( 13 ) 


Now let the square array with elements and i <_ j possess 
an inverse w’ ose elements are denoted by A 1 ^ ^ . Then this 
matrix may be used in (13) tc solve for the first M' sur- 
face densities, 
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, k = 1, M' 


PDEN 

(14) 


There are thus M-M' independent densities remaining and 
Eqn. (14) can be used to relate the dependent densities. 

The integration of the variational equations to 
obtair the partial s of the trajectory with respect to 
the independent surface densities requires that the fore 
ing function for the variational equations include both 
the direct and indirect effects of the independent ad- 
justed densities. If igjj is the surface density 
acceleration, then the required forcing function is 
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with - 7 ^— to bo obtained from Eqn. (14). 

It should be noted that GEODYN has the option of 
1 adjusting only a portion of the surface densities. This, 

in effect , means that there are additional constraint 
equations, but they are quite simple to incorporate. The 
constraints given by Eqn. (14) are still required to held 
with no modification whatsoever. Ordering the densities 
such that the unadjusted densities are last in the array, 
then Eqn. (15) is modified only to the extent that i has 
the range M’ + 1 to M-M 0 , with M q the number of unadjusted 
( densities. If there are more constraint equations than 

there arc densities to be adjusted, GEODYN will terminate 
upon reading the input deck with the appropriate e^or 
message . 





R.4 SOLAR . LUNAR. AND PLANETARY GRAVITATIONAL PERTURBATIONS 

SUNGRV 

The perturbations caused by a third body on a 
satellite orbit are treated by defining a function, 

R^, which is the third body disturbing potential. 

This potential takes on the following form: 





where 

m^ is the mass of the disturbing body. 

r d is the geocentric true of d isition 

vector to the disturbing bv 

S is equal to Jie cosine of the 
enclosed tie between r and r 

r is the geocentric true of date position 
vector of the satellite. 

G is the universal gravitational constant, 
and 

M is the mass of the Earth. 

The third body perturbations considered in GEOPYN 
are for the Sun the Moon, Venus, Mars, Jupiter, and 
Saturn. All areacompvted in subroutine SUNGRV by 
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where 
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These latter quantities, C and D as well as D 

are passed to subroutine VEYAL for the varia- 
tional equation calculations. VEVAL computes 
the matrix ^ whose i, j** 1 elements is given b/ 







( 3 ) 


This matrix is a fundamental part of the variational 
equations . 

8.5 SOLAR RADIATION MEASURE 

1 The force due to sola* radiation can have a 

significant effect on the orbits of satellites with 
a large area to mass ratio. The accelerations due 
to solar radiation pressure are formulated in the 

C 
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GEODYN system as 



Cl j 


where 

v is the eclipse factor, such that 

v=0 when the satellite is in the Earth's 
shadow 

v=l when the satellite is illuminated 
by the Sun 

is a factor depending on the reflective 
characteristics of the satellite, 

A g is the cross sectional area of the 
satellite; 

is the mass of the satellite, 

E s is the solar radiation pressure in the 
vicinity of the Earth, a_.d 

r s is the (geocentric) true of date unit vector 
pointing to the Sun. 

The unit vector r s is dete ined i s part of the 
luni- solar-planetary cphemeris computations 




The eclipse factor, v, is determined as follows: 


Compute 
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where r is the true of date position vector of the 
satellite. If D is positive, the satellite is always 
in sunlight. If D is negative, compute the vector P R 


P R - r - D r s . 


(3) 


This vector is perpendicular to r g . If its magnitude 
is less than an Earth radius, or rather if 






(4) 


the satellite is in shadow. 

The satellite is assumed to be specularly 
reflecting with reflectivity p s ; thus 


* 1 * p s 


( 5 ) 


When a radiation pressure coefficient is being 
determined; i.e., C R , the partials for the f matrix 
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in the variational equations computation must be 
computed. The i^ element of this column matrix is 
given by 
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f i * * v - P s r 

m s 1 
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These computations for the effects of solar 
radiation pressure are done in subroutine F. 


8.6 ATMOSPHERIC DRAG 


A satellite moving through an atmosphere ex- 
periences a drag force. The acceleration due to 
this force is given by 


DRAG 
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Cp is the satellite drag coefficient 


A. is the cross-sectional area of the satellite 
s 
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n? s .■ s the mass of the satellite, 


DRAG 


Pq is the density of the atmosphere at the 
satellite position, and 

v r is the velocity vector of the satellite 
relative to the atmosphere. 

Both A and C T . are treated as constants in GEODYN. 
s o 

Although A s depends somewhat on satellite attitude, the 
use of a mean cross-sectional area does not lead to 
significant errors for geodetically useful satellites. 

The factor varies slightly with satellite shape and 
atmospheric composition. However, for any geodetically 
useful satellite, it may be treated as a satellite 
dependent constant. 

The relative velocity vector, v r , is computed 
assuming that the atmosphere rotates with the Earth. 

The true of date components of this vector are then 


x r = i - e g y 


( 2 ) 


y r s y - e g x (3) 


i r - i (4) 

as is indicated from Section 3.4, the subsection on 
transformations between Earth* fixed and true of date 
systems. The quantities x, y, and z are of course the 
components of r, the satellite velocity vector in true 
of date coordinates. 
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* 

The drag accelerations arc computed in the 
GEODYN system by subroutine DRAG, with the atmospheric 
density Op being evaluated by subroutine D71, D650. in 
( addition, subroutine DRAG computes the direct partial s 

for the f matrix of the variational equations when the 
drag coefficient Cp is being determined. These partials 
are given by 
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When drag is present in an orbit determination 
run, the D r matrix in the variational equations must 
also be computed. This matrix, which contains the 
partial derivatives of the drag acceleration with 
respect to the Cartesian orbital elements, is con- 
structed in subroutine VEVAL. We have 
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where 
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is (x,y ,r.,x,y,x) ; i.e., x t spans r and F. 
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is the matrix containing the partial deriva- DENSTY 
tives of the atmospheric density with respect 
to x t and is partially computed in subroutine 
DENSTY (see section 8.7.4 on atmospheric 
density partial derivatives). Because the density 

is not a function of the satellite velocity, 

*Pn 


ar 


the required partial* are 
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8.7 ATMOSPHERIC DENSITY 

In the computation of drag, it is essential to li/i 

obtain models of the atmospheric density which will yield 
realistic perturbations due to drag. The GEODYN program uses 
the 1971 revised Jacchia Model which considers the densities 
between 90 km and 2500 km, and the 1965 Jacchia-Nicolet 
Model which gives densities between 120 km and 1000 km with an 
extrapolation formula for higher altitudes. 

The following discussion will cover primarily the assump- 
tions of the models and empirical formulae used in subroutine 
D71 and subroutine D650, The procedure for empirically 
evaluating the density tables which was developed by WOLF will 
also be included in the discussion. 

8.7.1 JACCHIA 1971 DENSITY MODEL 

The 1971 revised Jacchia model, as implemented in sub- 
routine D71, is based on Jacchia’ s 1971 report, ’’Revised Static 
Models of the Thermosphere and Exosphere with Empirical 
Temperature Profiles" (Reference 1) . The density computation 
from the exospheric temperature as well as from variations 
independent of temperature is based on density data appearing 
in that report. This data presented in Table 1 shows the density 
distribution at varying altitudes and exospheric temperatures. 

For further detailed development of these empirical 
formulae, the interested reader should consult the afore- 
mentioned report and Jacchia's 1970 report (Reference 2). 



8. 7. 1,1 The Assumptions of the Model 


The Jacchia 1971 model (J71) is based on empirically I 
determined formulae with some inherent simplifying assump- 
tions. Such an approach is taken primarily because the 
various processes occurring in different regions of the 
atmosphere are complex in nature. Moreover, at present, a 
thorough comprehension of such processes is lacking. The 
present J?i model is patterned after the J65a (Jacchia 1965a) 
model which was based upon previous assumptions by Nicclet 
(Reference 3) . 

In Nicolet's atmospheric model, temperature is con- 
sidered as the primary parameter with all other physical 
parareters such as density and pressure derivable from 
temperature. This approach was adopted by Jacchia in his 
J65a model. However, in the J71 model, there are variations 
modeled by Jacchia which are independent of temperature. 

They are the semi-annual variations, seasonal -latitudinal 
variations of the lower thermosphere, and seasonal - 
latitudinal variations of helium, all of which involve a 
time dependency. Although in J71 Jacchia mentions variations 
in hydrogen concentration, he does not attempt any quantita- 
tive evaluation of these variations. 

Composition 

The J71 model has assumed that the only constituents 
of the atmosphere are nitrogen, oxygen, argon, helium, and 
hydrogen. This composition is assumed to exist in a state 
of mixing at heights below 100 km and in a diffusion state 
at higher altitudes. A further assumption on the composition 
of the atmosphere is that any variation in the mean molecular 
mass, M, in the mixing region is the result of oxygen dissocia- 
tion only. The variation in M has been described by an empiri- 
cal profile for heights ranging from 90 km to 100 km. 
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It is also believed that gravitational separation 
J for helium exists at lower height than for the other compo- 

nents. To avoid the cumbersome ordeal of modeling a separate 
homopause for helium, Jacchia has modified the concentration 
at sea-level by a certain amount such that at altitudes 
^ where helium becomes a substantial constituent, the modeled 

densities will correspond to the observed densities. 

Although this will yield a higher helium density below 
100 km, the contribution of helium to the overall density 
will be negligible below this height. 

\ 

Hydrogen does not become part of the density model 
until a height of 500 km. At this altitude, hydrogen is 
assumed to be in the diffusion equilibium state. 

Temperature 

The temperature above the thermopause is referred to 
^ as the exospheric temperature . Although this temperature 

is independent of height, it is subject to solar activity, 
geomagnetic activity, and diurnal and other variations. 

The J71 model assumes constant boundary conditions of 90 km 
with a constant thermodynamic temperature of 185° K at this 
height. From numerous atmospheric conditions it is suggested 
that the atmospheric conditions at 90 km do indeed vary 
nominally, and thus, this assumption may be reasonably 
acceptable (Reference 4) . Profiles of the thermodynamic 
temperature show that it increases with height and reaches 
an inflection point at 125 km. Above this altitude, this 
temperature asymptotically attains the value of the 
exospheric temperature. An analytic model of the atmospheric 
densities by Roberts (Reference 4) has been constructed based 
on modifications to Jacchia' s 1970 temperature profile between 

( 
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90 km and 125 km The J71 model assumes that the basic 
shape of the temperature profiles remain unchanged during 
atmospheric heating due to geomagnetic storms. In all 
liklihood, the profiles at low altitudes become distorted 
to yield higher temperatures during such occurrences. 

Since the J71 model assumes the atmosphere to be in 
static equilibrium, for any sudden changes in solar activity 
or in geophysical conditions, which are characteristically 
dynamic, the model will generally be unable to properly 
represent the variations in both temperature and density 
due to this invalid assumption. The priority has "been 
given to the best representation of density. 

8 . 7 . 1 . 2 Variations in the Thermosphere and Exosphere 

Several types of variations occurring in the different i 

regions of the atmosphere are incorporated in the J71 model. ! 

1 

These variations are: solar activity variations, diurnal 

variations, geomagnetic activity variations, semi-annual 
variation, seasonal-latitudinal variations of the lower 

| 

thermosphere, and seasonal-latitudinal variations of helium. j 

Still another variation which is not quantitatively evaluated 
by J71 is the rapid density fluctuations believed to be j 

associated with gravity waves (Reference 1) . Each of the I 

above variations may be modeled empirically from observable j 

data. However, because a static model is used, the various ! 

predictions will exhibit different degrees of accuracy for j 

each variation. It is fundamental, however, to note that j 

unless the characteristic time for which these variations j 

occur is much longer than that for the processes of diffusion, 
conduction, and convection to occur, the predictions will be 
poor (Reference 1) . 
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Solar Activity 

The variations in the thermosphere and exosphere as 
a result of solar activity are of a dual nature. One type 
of variation is a slow variation which prevails over an 
11-year period as the average solar flux strength varies 
during the solar cycle. The other type is a rapid day-to-day 
variation due to the actively changing solar regions which 
appear and disappear as the sun rotates and as sunspots are 
formed. 


To observe such activities, the 10.7 cm solar flux 
line is commonly used as an index of solar activity. The 
National Research Council in Ottawa has made daily measure- 
ments on this flux line since 1947. These values appear 
monthly in the "Solar Geophysical Data (Prompt Reports)" by 
the National Oceanic and Atmospheric Administration and the 
Environmental Data Service in Boulder, Colorado (U.S. 
Department of Commerce) . 

A linear relationship exists between the average 10.7 
cm flux and the average nighttime minimum global exospheric 
temperature (Jacchia 1971) and may be expressed as: 


= 379° + 3.24® F 10>? (’Kelvin) (1) 

where 

T ro = is the average nighttime minimum global 

exospheric temperature averaged over three 
solar rotations (81 days) . 
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F 10.7 is the average 10.7 cm flux strength over 

three solar rotations and measured ir. uni\.r 
of 10 watts m’^ (cycle/sec) handwic* . . 

Equation (1) expresses the relationship with solar Tlux when 
the planetary geomagnetic index, Kp is 2 eror, i.e., for no 
geomagnetic disturbances. 

The nighttime minimum of the global exospheric tempera- 
ture for a given day (Reference 1) is 


T c + 1,3 ° ^ F 10.7 ' ^10.7^ W 

where 

F^q y is the daily value of the 10.7 cm solar flux 

in the same units as F 1Q y for one day earlier, 
since there is a one day lag of the temperature 
variation response t.o the solar flux (Roemer 
1968) . 

Thus, Equation (2) models a daily temperature variation 
about the average nighttime minimum global temperature as 
determined in Equation 1. 

Diurnal Variations 

Computations from drag measurements, have indicated 
that the atmospheric density distribution varies from day 
to night. The densities axe at a peak at 2 P.M. local 
solar time (LST) approximately at the latitude of the sub- 
solar point, and a minimum at 3 A.M. (LST) approximately 
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of the same latitude in the opposite hemisphere. The 
diurnal variation of density at any poin* is ouoject to 
seasonal changes. By empirical relationships, this varia- 
tion may be described in terms of the temperature. Again, 
because a static model is used, the accuracy of this 
temperature is open to question. 

At a particular hour and geographic location, the 
temperature, T^, can be expressed in terms of the actual 
global nighttime minimum temperature, T , for the given 

v 

day (Reference 1). Thus, we may write 
T £ » T c (i + R sin m 6) ^1 + R 

where 

R = 0.3 
m * 2.2 
n = 3.0 

t » H + 0 *• p sin(H+y) for ( 

P» * -37° (lag of the temperature maximum with the 
uppermost point of the sun.) 
p « + 6° (introduces an asymmetry in the temperature 
curve . ) 

y " +43° (determines the location of the asymmetry 
in the temperature curve.) 

n - j abs (♦’-J,) 
a - £ abs (*'♦«,) 


m . m„ \ 

cos n - sm e t V 

= COS - 1 

1 + R Sin m e 2 f 


(3) 
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<J> ' * geographic (geocentric) latitude D"1 

6.= declination of the sun 

H - hour angle of the sun 

(when the point considered, the sun, and the 
earth's axis are all coplanar, H=0. The hour 
angle is measured westward 0° to 360°.) 

The parameter R determines the relative amplitude 
of the temperature variation. Jacchia and his associates 
have undertaken investigations of R which reveal indications 
of its variation in tine and with altitude. After consult- 
ing 1969-1970 data, Jacchia presently has abandoned any 
attempt at any definitive conclusions about the variations 
of R with time or with solar activity (Reference 1) . In- 
stead, he believes this evidence to be the result of inherent 
limitations of the static atmospheric representation. Con- 
sequently, in the J71 model, a constant value of R-0.3 is 
maintained. 

Geomagnetic Activity 

Precise effects of geomagnetic activity cannot be 
obtained by present measurements from satellite drag, since 
suca techniques can only show averaged values of densities. 

It is also realized that the consequences of a geomagnetic 
disturbance in view of the atmospheric temperatures and 
densities over the global regions are of a complex nature. 
However, when such disturbances occur, there are indications 
of increases in temperature and density in the thermosphere 
above the aurora belt. By the t*'i« this atmospheric dis- 
turbance reaches the equatorial . • o.is, a period of roughly 
7 hours, the effects are damped . .t considerably. (ReTerence 1) . 
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A static model description of tempera cure and density 
cannot be viewed accurately since the propagation time of 
the geomagnetic storms is relatively - v ort. Therefore, any 
empirical formulae used to compute the effoetr cn the para- 
meters yield only a vague picture. 

Jacchia et al_ (1967) have related the exospheric 
temperature to the 3 -hourly planetary geomagnetic index K p . 
quantity K p is used as a measure of a three -hour variation in 
the earth’s magnetic field. The response of the temperature 
change to geomagnetic storms lags the variation in K p by 
about 6.7 hours. In the following equation (Reference 1) the 
correction to the exospheric temperature due to geomagnetic 
activity is 


AT. - 28° Kp + 0.03° exp (K p ) 


( 4 ) 


for heights above 200 km. 

Although this K p in equation (4) is a three -hour 
planetary geomagnetic index, in subroutine DENSTY an 
averaged K p over a 24 -hour period is used to minimize the 
amount of *nput data to GEODYN. The loss of accuracy in 
using the daily average of K p is minimized, since the above 
equation is for a smoothed effect of the variations derived 
from satellite data. 

Below 200 km, density predictions from equation (4) 
prove to be too low. Better results are obtained if the 
variations were represented as a two-step hybrid formula 
in which a correction to the density and to the temperature 
is made. Thus, in J71 the following hybrid formula 
(Reference 1) is given as 


The 
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(a) P 4 = Alog 10 P - 0.012 K ? + 1.2 x 10' 5 exp (K p ) 

(b) AT W - 14° Kp ♦ 0.02 e exp (X p ) ( 5 ) 


where Alog 1Q p is the decimal logarithm correction to the 
density p. 

The values of a three -hour K p index are available 
along with the daily solar flux data in the publication 
'.'Solar Geophysical Data" by the N.O.A.A. and E.D.S., 
Boulder, Colorado (Department of Commerce). 


In computing the exospheric temperature, accurate 
daily values for both the solar and geomagnetic flux must 
be used. These values are stored in the subroutines FLUXM 
and FLUXS of GEODYN, and they are updated as new informa- 
tion is received. This information may be updated (sub- 
routine ADFLUX) using the appropriate GEODXN Input Cards. 

The user should be aware of the fact that these tables are 
expanded as new information is made available (Reference 3) . 


FLUXM 

FLUXS 

ADFLUX 


At the beginning of each run, a file is generated for JANTHG 
each satellite arc which contains the required flux data 
for the time period indicated. Subroutine JANTHG sets up 
the flux tables as well as averaging the daily values of 
solar flux over three solar rotation periods. This enables 
the releasing of vast computer storage required for daily 
flux values over 14 years. The selected data is stored in 
common block FLXBLK. 


A midpoint point average is used to compute the six 
solar rotation flux values 
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Semiannual Variation 


The semiannual variation at present is least under- 
stood of the atmospheric variations. In past models, J65, 
attempts at empirically relating the temperature to this 
variation seemed appropriate in the range of heights, 

250 to 6S0 km, for which data was available. However, with 
tat availability of new data for a wider range of altitudes, 
laig r discrepancies in the densities appeared. After close 
scru'lny, Jacchia in 1971 (Reference 1) found that the 
’ litudeci the semiannual density does net appear to be 
connected with the solar activity. It does, however, show 
a strong depend r„ce on height and a variation from year to 
year. Drag analyses from the Explorer 32 satellite have 
revealed that a primary minimum in July and a primary maximum 
in October occur for the average density variation (Reference 
1 ). 


Jacchia in J71 expresses the semiannual ariation 
as a product function (Reference 1) in which 


P 2 » Alog 10 o - f(z)g(t) 


( 6 ) 


where f(z) is the relationship between the amplitude, i.e., 
the difference between the primary maximum and minimum, and 
the height, z, and where g(t) is the average density varia- 
tion as a function of time for the amplitude normalized to 1. 
The two expressions for f(z) and g(t) which yield the best 
fit to the data are 


8.7-11 


371 


f(z) = (5.876 x 10' 7 z 2,331 ♦ 0.06528) exp (-2.868 x 10" 3 z) 

(7) 


for i in kilometers; 

g(t) * 0.02835 + 0.381711 + 0.4671 sin (2Ih + 4.1370) 
sin (4n + 4.259) 

( 8 ) 

where 


t - 4 + 0 . 09544 1 [ 0 . 5 + 0.5sin (2ir$ + 6.035)] 1 * 650 -O.sj 

# = (t - 36204) /365 . 2422 

t * time expressed in Modified Julian Days 
(M.J.D. ■ Julian Day minus 2 400 000.5). 

M.J.D. c 36204 is for January 1, 1958. 

The term * is the phase of the semiannual variation which 
is the number of days elapsed since January 1. divided by 

the number of days for the troplr.#' 1 year. 

Seasonal -Latitudinal Variations of the Lower 
Thermosphere 

In the lower thermosphere, from about 90 km to 120 km, 
there are variations occurring in temperature and density 
depending on the latitude and the season. Only the density 
variation is considered because any temperature variation 
proves to be too tedious to handle. Between the heights 
from 90 km to 100 km, there is a rapid increase in the 
amplitude of this variation in density with a maximum ampli- 
tude occurring between 105 and 120 km (Reference 1). Above 



D71 


120 km there is no data on which to base prediction? of 
the seasonal-latitudinal variations. This variation 
appears to decrease in amplitude to the point where 
negligible fluctuations exist at ISO km. Therefore, in 
DENSTY, the seasonal-latitudinal variations are neglected 
at heights above 160 km. 

Jacchia in J 7 1 fits the seasonal variations tc an 
empirical correction to the decimal logarithm of the 
density (Reference 1) as follows: 

P, = Alog in p * S P sin 1 +' (9) 

* 1 U’l 

where 

i»« * geographic latitude 

S = 0014 (z-90) exp [ -0 . 0013 (z - 90) 2 ) 

z = height in kilometers 

P = si*i + 1.72) 

$ « phase as given in equation (8), 
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aeasonal-Latitudinal Variations of Helium 

Helium in the atmosphere has been observed to n»i2 rat - e P/1 
towards the winter pole. The phenomenon of this seasonal 
shift in the helium concentration in the upper atmosphere 
is not yet understood. It therefore becomes necessary to 
perform an empirical fit to drag data from which this 
seasonal variation is derived. The expression which is 
used in J71 (Reference 1) to describe the helium variation 


I | r r / ^ $ 

Q 2 = -log 10 n (He) = 0.65 j — J in ' ^ - - - 


2 | 6 »| 


- sin'' — (10) 

4 J 


where 


n(Hi) = number density of helim (number of particles/cm ) 
A© = declination of the sun 
= geographic latitude 

e = obliquity of the ecliptic (e = 23.44°) 

The variation of the helium density in subroutine 
DENSTY is not considered for heights below 500 km. It is 
also neglected for latitudes whose absolute value is less 
than 15° between the range of heights from 500 km to 800 km. 


The correction to the density due to the seasonal lati- 
tudinal variations of helium is then 


Ap n * 10 


]og, 0 n( H e)f Alog 10 n(He) ] 

1U 10 1U -1C gm / cm 3 


where 


C is the molecular mass of Helium divided by Avogadro’s 
Number. 
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8, 7, 1.3 Polynomial Fit of Density Tables 


The data which appears in Table 1 shows the variation D71 
of density with altitude and exospheric temperature which 
is reproduced from Jacchia’s 1971 report (Reference 1). 

From heights of 90 km to 100 km, the density values weie 
obtained by numerically integrating the barometric equa- 
tions. The diffusion equation was numerically integrated > 

to obtain values of the density on the altitude range, i 

100 km. <Z< 2500 km. In both cases, an empirical tempera- j 

ture profile was used for each exospheric temperature. j 

In the GEODYN subroutine DENSTY, the atmospheric 
density is computed based on the data from Table 1 after < 

appropriate corrections are applied to the exospheric tempera- 
ture. The tabulated densities have been fitted (by WOLF) to 
various degree polynomials of the form: 


P 1 - ^lo’OT ■ £ 

i j 

where 

Ppj is the density in g/cm 3 

T is the exospheric temperature, 

h is the spheroidal height (altitude) , and 

•jj is a set of appropriate coefficients for 
tho density tables. 
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third degree fit, the coefficient# for the selected polynomials 
fnr the total dens’. V *.'* shown in Table 2. In Table 3, coeffi- 
cients of polynomials for the helium number density uic pic 
sented. 

The commuted densities from the fitted polynomials D71 
show a reasonable percentage error from the densities given 
in Table 1. For each of the regions and temperature ranges, 
the maximum errors are given in Table 4. The largest error 
of 12% occurs in the reg5 n n between 500 - 1000 km in the 
temperature range of 500° - 800°K. In the region of 
1000- 2500 km with temperatures between 800° • 1900°K, a 
fourth degree fit to the temperature yields a maximum 
error of 11.0% in the densities. 

The helium number density fits are also given in 
Table 4. As one can see, the values of the number density 
are quite satisfactorily fitted by the polynomials. The 
maximum error in the whole range of heights and temperatures 
is only 2.8%. 

Overall, these fits could be improved by either using 
higiier degree polynomials or possibly other functions, or 
by further sub-dividing the density table. However, these 
maximum ei rors appear to be tolerable since they are con- 
sidered to be within the range of accuracy of the model 
presently used. Above 2500 km, the density was found to 
be negligibly small, and therefore, was set to zero. 
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Summary of Log Densities 
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TABLE 2. 

DENSITY POLYNOMIAL COEFFICIENTS 
(For Decimal Logarithm of Density) 



T° 

T 1 

T 2 

T 3 

90-200KM 

h° 

4.22085 

0 . 98393E-2 

- . 64952E- 5 

0 . 1471SE-8 

h 1 

-0.20134 

- . 23412E-3 

0 . 1S337E-6 

- . 34675E-10 

h 2 

0.78592E-3 

0. 16966E- 5 

-.S1060E-8 

0 . 25007E-12 

h 3 

- . 12087E-5 

34360E-8 

0.22457E-11 

51069E-15 


200-500KM for 500° -800°K 

h° - . 12838E+2 0.40709E-2 
h 1 0. 82282E-1 - . 3121 SE-3 
h 2 - . 689S1E- 3 0. 24402L- 5 
h 5 0.11263E-5 - . 41807E- 8 


0.97074E-5 -.10643E-7 

0. 26S43E-6 -.SS193E-10 

-.270S8E-8 0.99003E-12 

0.50617E-11 - .20484E-14 


200- SOOkM 



for 800°-1900“K 
-3.4595 - . 15000E- 3 

- . 28395E- 1 0 17760E-6 

0. S5998E-S 0.77461E-7 

0.39434E-8 -.76435E-10 


- .62640E-.6 0. 24612E-9 

0 61398E-8 - . 23362E- 11 

- . 59492E-1C 0.14921E-13 

0. 58333E-13 -.14S9SE-16 


500-100KM for 500 # -8Q0°K 

h° - . 77659E+2 0.167271 -.56S70E-4 -.50424E-7 

h 1 0.50638 - .98936E-3 0.749S2F.-6 -.55178E-10 

h 2 - . 38935E-3 0.12973E-S -.19776E-8 0*1419111-12 

h 3 0. 1E962E-6 -,$4949E-9 4>.4670*E-12 -.7I8tt6E-XA 
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TABLE 2 

DENSITY POLYNOMIAL COEFFICIENTS 


m I 


500-1000KM for 800°-1900*K 
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0.21370E-16 
- .41504E-20 



TABLE 3 


HELIUM DENSITY POLYNOMIAL COEFFICIENTS 
(DECIMAL LOG OF HELIUM NUMBER DENSITY) 


T 1 


500-1000KM for S00°-80C°K 

h° 9.3712 52634E-2 

h 1 - . 13141E-1 0. 31218E-4 

h 2 0. 26071E- 5 - . 7S730E-8 

h 3 - . 521S6E-9 0 . I9056E-11 


0. 52983E-5 
- . 32593E-7 
0. 83058E-11 
- . 26578E-I4 


-. 20471E-6 
Q.12S73E-10 
- .40669E-14 
0.12S3SE-17 


SOO-IOOOKM 

for 800°01900 

i°K 

h° 

8.3914 

- . 16433E-2 

h 1 

- . 69049E- 2 

0.84138E-5 

h 2 

0. 1O510E-5 

- .12663E-8 

h 5 

- . 12222E-9 

0. 14745E-12 


0. 78032E-6 
- .44577E-8 
0. 71134E-12 
- .97658E-16 


-.14323E-9 
0.85627E-12 
- . 14180E-15 
0. 21458E-19 


1000-2500KM for 500° -800 
h° 9.1045 

h 1 - . 12259E-1 

h 2 0. 15893E-5 

h 3 - . I1829E-9 


*K 

-4. 3410E-2 
0. 27951E- 4 
- . 35863E-8 
0. 26138E-12 


0.40292E-5 
. - . 27972E-7 
0. 3S476F.-H 
- . 25227E-15 


-.14S22E-8 
9.10371E-10 
- .1298SE-14 
0.89714E-19 


100-2500KM for 800*-1900 # K 

h° 8.6120 -. 25363E-2 

h 1 - , 84847E-2 0.14084E-4 

h 2 0.11543U-5 -19884E-8 

h 3 -.94521F.-10 0. 17387B-12 


. . ' 


0 . 18979E- 5 - . 73696E-9 0.1138$E-T2 

- .11386B-7 0.44871E-J> 

0 * 16635E-11 - . 0 } 

-.15368C 




TABLE 4 


f* 




PERCENTAGE ERROR OP POLYNOMIAL FITS 
TO THE DENSITIES 


( 

Height 

Range 

(KM) 

Temperature 

Range 

(°K) 

Maximum Percent Error 
Total Density Helium Density 


90-200 

500-1900 

11.0% 

- 

( 

200-500 

500-800 

11.6% 

- 


200 500 

800-1900 

5.13% 

- 

‘ 

500-1000 

500-800 

12.0% 

0.441 

f 

500-1000 

800-1900 

8.85% 

2.8% 


1000-2500 

500-800 

4.1% 

1.01 


100-2500 

800-1900 

11.0% 

1.251 


( 




■ +* 


8 . 7. 1.4 The Density Computation 

When all of the terms contributing tc the atmosphere 97 ^ 
density are combined 


3 I ^l + ^2 + ^3 + ^4 ^1 ^2 

P D “ 10^ 10 1 c 0 * * 10 A (10 L - 1) C 


(14) 


where 


p n * the atmospheric density in Kg/m' 


Pj is given by equation (12), 
is given by equation (6), 


Pj is given by equation (9) , 

P4 is given by equation (5a), 

Qj is given by equation (13), 

Q2 is given by equation (10), and 

C is the molecular mass of Helium divided by 
Avogadro’s Number « 0.6646(10“^) 
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8. 7. 1.5 Density Partial Derivative:; D71 

- VUVA 

JD- 

$ I ri addition to the density, GHODYN also requires the 

partial derivatives of the density with respect to the 
Cartesian position coordinates. These partial-s are used 
in computing the drag contributions to the variational 
equations. 

The spatial parti 1 derivatives of the atmospheric 
density are 


( 



3p D _ 

3P0 

3* 

+ 

3p D M 

_ <f 

ifo 

3h 

3r 

3$ 

3r 

3* 

3h 

3r 


(i) 


where 


h - spheroid height of the satellite 
$ - sub-satellite latitude 
X - sub-satellite longitude 

r - true of date position vector of the satellite 


Variations in atmospheric density are primarily 
due to changes in height. Therefore, only height varia- 
tions are computed by CEODYN’ and 


3* 

3X 


( 2 ) 


0 



and consequently 


J71 


where 


where 


3Pq 3 Py Sh 

3r 3h 8r 


(3) 


ah 

Sr 


is presented along with the spheroid height 
computation in Section 5.1. 


The density is given (Section 8.7.4) by 


P D * 10- 


10 


P l +P 2 +P 3 +P 4 « 


+ 10 A (10 * - 1 ) c 


(4) 


Pj-j - density in Kg/m' 


n 


m 


p i - E h(i ' 1) E tU ' 1) 


i=l 


(5) 


P 2 “ g(t) [5.876(10“ 7 ) h 2 * 3jl + 0.06328] exp [- 2. i68 (10” 3 )h] 


♦"2 2 
P- *0.014 (h-90) P sinV exp t-0.0013(h-90r] (7) 

3 . in 


-5. 


P. = 0.012 K + 1.2(10 a ) oxp (K ) 

P P 


( 8 ) 
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D71 


( 


n 


m 


Q i = E E b ij t °'” - io sio "( |ic ) 


(9) 


i = l 


j=l 


Q 2 = ^log 10 n(Mc) 


( 10 ) 


C = the molecular mass of Helium divided by Avogadro's 
Number. 

h = height in Km. 


- polynomial coefficients used to fit. the density 
table. 


**ij “ polynomial coefficients used to fit the Helium 
number density table. 


( 


All other terms are defined in Section S.7.1,4 and need no 
further clarification at this point since they are constants 
in the partial derivative equations. 


Defining two basic derivative .formulae , 




1 ,»(x) . ,»(*) i uW . 

dx 


dx 


( 11 ) 


* f 

i 


: I 
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D71 


t. lo u(x) = i. c u(x)*nlO 


lulO c 


u(x) AnlO 


du(x) 


c io u ^ x ^j>nlO 


du(x) 


( 12 ) 


And it follows that 


3 p i +p 2 +p 3 +p 4 p i +p 2 +P 3 +P 4 3 

— 10 1 2 3 4 = 10 1 2 3 4 InlO — (P.+P-+P-+P.) 

Oh Oh 1 1 3 


(13) 


0 Q, Q, OQ, 

— 10 1 « 10 1 JlnlO — - 

Oh Oh 


(14) 


Differentiating the components of (13) and (14) 


*s: 


> it 


(15) 


= g(t)|:.B76(10“ 7 ) (2.331) h 1,331 exp [-2. 868(10"' 
♦ [S. 876(1 O’ 7 ) h 2,331 + 0.06328] (-2.868) (10' 3 ) 


exp [-2.86S(10 J ) h) > (16) 
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M**" *( 


0.014 P 


sin 2 *' exp [-0. 0013(h-90) 2 ] 


*-> t x 


3P. 




ah 




2 


1 + 2 (h- 90) “ (-0.0013) 


! 


(17) 


3P, 

ah 


*= c 


(18) 


3Q 

ah 


m 


i - 2 U*« h(i ' 2) L a ij T °' 1) 

j-1 


i=2 


(19) 


c 


3 

The resulting partials are in the units of (Kg/m )/Km 

-3 


and must therefore be multiplied by 10 


apjj a 
ah ah 


^ > l + ^ > 2 + ^ > i + ^ > 4 ^2 w xi 

10 1 C * q + (10 L - 1) C — 10 1 


8 Q) 
1( 
ah 


VEYAL 

( 20 ) 


The units of (20) are then 


(Kg/m 4 ) . 
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8 . 7 . 2 JACCI1I A 1965 Density Model 

D65C . 

The J tech in 1965 Density Mode], as implemented in sub- 
routine D650 is based or Jacchia's J965 Teport, "Static 
Diffusion Models of the Upper Atmosphere with Empirical Tempera- 
ture Profiles" (Reference 12). The formulae for computing the 
exospheric temperature have in some cases been modified according 
to Jacchia's later papers. The density computation from the 
exospheric temperature is based on density data provided in that 
report, reproduced herein as Table 5, which presents density 
distribution versus altitude and exospheric temperature . 

The reader who is interested in the development of these 
empirical formulas and the reasoning behind them should consult the 
above mentioned report and Jacchia's later papers. For the con- 
venience of this interested reader, the references 13 for this sec- 
tion from a reasonable comprehensive bibliography. 
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8. 7. 2.1 The Assumptions of the Model 

• The lacchin-Nicolct model is based on certain simplifying 
assumptions and on empirically determined formulae. This is pri- 
marily due to the complexity and varied nature oT the processes 
occurring in different regions of the atmosphere and the general lack 
of anything rescnbl ing a complete understanding of the fundamental 
mechanisms involved. The actual derivation of the model is based 
upon assumptions first proposed by Nicolct (see Reference 14); Jacchia 
selected the Nicolet approach to generate a model suitable for 
satellite dynamics. 

The model of the atmosphere proposed by Nicolet considers that 
the fundamental parameter is the temperature. Other physical parametc 
such as the pressure and density were derived from the temperature. 
Thus the first concern is the temperature variation in the atmosphere. 

This temperature variation is controlled by the following 
conditions: 


1, Above the thermopause, the temperature of the atmosphere 
does not vary with altitude. The thermopausc varies with 

m 

solar activity (and the time of day) , ranging 
between about 220 km to 400 km. The 
temperature above the thermopausc is called D650 

the exospheric temperature and is directly 
responsive to solar effects. 
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D650 

2. At an altitude of 120 kin , the temperature, 

density, and atmospheric conditions arc inde- 
pendent of time. This is an obvious simpli- 
fication. Howcvc. , the variations of these 
parameters above 120 km arc considerably 
larger than those occurring at 120 km, and, 
considering the other assumptions, this 
assumption represents a reasonably good 
approximation. ‘ . 


3. The atmosphere is assumed to be-in static 
equilibrium. With the large day-to-night 
temperature variations, having a period oi the 
sane order of magnitude as the conduction time 
in the lover thermosphere, and vrith the oc- 
casional occurrence of severe magnetic storms 
which give vise to fairly rapid and large 
temperature variations the validity of this 
assumption is open to question. The best 
argument for this assumption is its relative 
simplicity. It should be anticipated, however, 
that in times of rapid change of the solar or 
geophysical parameters the predictions of this 
model will be in error due to the invalidity 
of this assumption. 


The atmosphere is considered to be in diffusive 
equilibrium above 120 km; that is, the density distribu 
tions of each atmospheric constituent with height arc 


7 . 30 
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governed independently by gravity and temperature. The 
governing equations arc the hydrostatic law, relating 
the pressure variation with height to the acceleration 
of gravity, and the perfect gas low, which relates the 
pressure, density and temperature. 

With this approach, Nicolct showed that above 
250 km the observed density profiles were reproduced 
satisfactorily if the (exospheric) temperature was as- 
sumed to be a different constant'. He also indicated 
that the problem of representing the density between 
120 km and the thermopause was largely a problem of de- 
ducing the vertical distribution of temperature. 


The contribution of Jacchia to the so-called 
Jacchia-N'icolct model is largely the development of 


empiric h 1 formulas to compute both the exospheric 


temperature and vertical temperature dir*i. uution as a 
function of exospheric temperature. These formulae are 
based on satellite observations coupled with physical 
reasoning. In addition, Jacchia has updated the boundary 
conditions of Nicolct. Thus in effect Jacchia has pro- 
vided all but the basic assumptions behind the model. 


The fundamental parameter of the model is therefore 
the exospheric temperature. This temperature, together 
with the boundary conditions, implies a particular vertical 
temperature profile. There three items - exospheric 
temperatur. , boundary conditions, and temperature profile - 
define the density at any altitude over 120 km through the 
diffusive equilibrium equation. 
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Figure 5, which was taken fror.. Reference 14, shows 
a comparison of density anu exospheric temperatures de- 
rived from observations of Explorer I satellite with 
solar and geomagnetic parameters . Note the correspondence 
between the exospheric temperature and the density. 


U6S0 


8.7.2. 2 The Exospheric Temperature Computations 

To calculate the fundamental parameter 3 the exo- 
spheric temperature, Jacchia considered four factors which 
could cause variations: . 

• «- 

• « 

X, Solar activity variation . 

• » 

2. Semi-annual variation 


Diurnal variation 


&E?RODUCIBILITY OF THE 
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4. Geomagnetic activity variation 

Each of these variations was determined to be related 
to one or more observable parameters (see Figure 1). 

The fc.ji.vcn empirical formulae are based on these parameters. 

Solar Activity 


There arc many indices of solar activity but the one 
whose variations most closely parrllcl those ol ntwj^crir* ;;- r 
density is the 10.7 cm. (2S00 Me.) boler'-tiwi f: 

.... intensity. of this line has been **nshi*^ 

'• ■ 5$4T, ; by . tlvO' : XttUo»a3; 'Research . 

Vssis„*Tl»c -vali*cs *f this 10.7 cn. 


v, ' 


F‘. J 


' i-.v-wss s-st j-: 
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Most of the time solar activity is much more intense 
in one solar hemisphere than the other so that the flux line 
appears to vary with the rotation period of the sun, 27 
days. This periodicity frequently persists for a year or 
•longer. In addition, there is a variation in the average 
flux strength Kith a period of about il years which is 
related to the solar cycle. ' 

• «. 

From satellite drag data a linear relation between 
the average 10.7 era. flux and the average global nighttime 
minimum exospheric temperature has been obtained (Reference 12) 

m * 

and is expressed as ' 

C 


monthly in the "Solar-Gcophys icn) Data Reports*' of the 
Environmental Science Services Administration in Boulder, 
Colorado (U.S. Department of Conmercc) . 


T # - *S7* ♦ 3. 60“ r i0> , 


v\ 


where 


10.7 
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is the average 1C,7 cm. flux strength over 
2 or 3 solar rotations measttred iit units 

of l(f 22 \vatts/n 2 /cyclc/sce, bandwidth, .. 

• v ■: •?;■••• - -v.v ■ • 

-/“■v ''l :• : . y , 4‘- “V; * * <v - - * - ./ 

is - the */; 


This formula gives tlie 1 ^ - 


*■ , * ;■ •<* 
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TIjc variation within one solar rotation is ex- 
pressed (Reference 12) by 


T 1 = T 


+ 1.8° (P 


10.7 *10.7 


( 2 ) 


where 


10.7 


is the mean cf the 1C. 7 cm solar flux 
for a given day in the same units as 

F 10 .7- and 

• » 
is the global nighttime minimum for the 
same day. 


This formula accounts (approximately) for the day to day 
temperature variation superimposed on the average global 
nighttime minimum temperature determined by the previous 
formula. 

There is some indication that, the coefficient 1.8° 

actually varies from sunspot maximum to sunspot minimum. 

The indicated range of variation is from about 2.4 P down 

to 1.5°. - f 
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Semi-Annual Variation 

The semi -annual variation is €Hc least understood, 
of the several types > 

sphere . Yearly , th‘0 -'tOO;'- \m, 
reaches a deep mi ni^sr'fld Ju t f issli 
maximum in Octobci>>lp^f§^^^ ..-;yv 

January,.; 
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(Reference 15) found that the observed density variations D650 
could be c. plained by temperature variations in the therno- 
pause, and arc roughly proportional to the 10.7 cm flux 
line. It /.as been noted that the height of the iono- 
spheric f 2 layer shows a semi-annual variation almost 
exactly in phase with the observed density variations. 

Another suggestion by P.S. Johnson (Reference Id) concern- 
ing the cause of the semi-annual variation, involves 
convective transfer at ionospheric levels from the 
summer pole to the northern pole. This, as yet, does 
not seem to account correctly for all the details 
of this variation. The semi-annual variation is not 
as stable a feature as the diurnal variation. Jacchia 
(Reference 12) accounted for this feature in 1965 but has, 
with the recent information of drag data from six satel- 
lites, updatedhis empirical formula (Reference 6) as 
follows : 


T 0 * T 0 + 2 ‘ 41 * I r 10t7 l«‘349+0.206 sin(2irt+226.5 c )] 

( 3 ) 


sin (4ttt+247.6°) 


where 





d 


day of the year counted from January 1. 


D6 


Y ® the tropical year in day.",. 

Tq - global night tine minimum temperature for 
that day corrected for semi-annual varia- 
tion. 

Jacchia, Slovcv, and Campbell (Reference 17) have more 
clearly defined this variation. As expected, the re- 
lationship between the temperature and the 10.7 cm flux 

* 

line cannot be considered accurate. It was concluded 
that the observed density variations are the result of 
temperature variations at essentially the same level as 
in the case of the solar effect. However, a variable 
altitude shows that the semi-annual variation affects 
the whole .atmosphere in the sane manner, irrespective 
of latitude. 
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Diurnal Variation 


D650 


The most regular of the variations is the diurnal 
variation. One can picture the density distribution as 
an atmospheric bulge with its peak 30° cast of the sub- 
solar point, degrading nearly symmetrically on all sides, 
but a little steeper on the morning side. The density 
peaks at 2 P.M. local solar tine and the minimum occurs 
at 4 A.M. The ratio of the maximum temperature 
at the center of the bulge to the minimum in the opposite 
hemisphere remains constant throughout the solar cycle; 
the ratio is 1.28 in Jacchia's model, atmosphere . ‘ The 
cause of the heating is in dispute. Some investigators 
believe it is due entirely to extreme ultra-violet (EUV) 
radiations; others, to ion drift; and still others, to a 
combination of the tv:o. 


The temperature, T, at a given hour and geographic 
location, can be computed in terms of the correct global 
nighttime minimum temperature for that day, Tq, using 
the following formula which approximates a mathematical 
description of the atmospheric bulge (Reference 12): 


T q (1+R sin 0)1 1 + - 


RCcos^n-sin 1 ^) t \ 

cos - ] 

1+R siii 0 2 f 
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where 


D053 


R *■ 0.26 


n = m = 2 . S 


t = 11 + B + p sin (H+y) (-u<t<n) 


B = -45 e 


12 ( 


V « 45’ 


AbSlft-' « e )/2] 


A 

v 


l f ^ f N 

i / <. j 


-geographic latitude * 

REPRODUCIBILITY OP DIE 

declination of the sun ORIGINAL PAGE IS POOR 

H c hour angle of the sun 

(H « 0 occurs when the point considered, 
the sun, and the earth’s axis are coplanar, 
II is measured westward 0° to’ 360°) 


Based on satellite information, dncchia (Reference IS) assumes 
9 maximum day tesaperaturo 28 * higher than the corres- 
ponding nighttime minimum. The variation is represented 
by R in the above equation. However, further investi- 
gation by Jacchia, Slou'ey, and Campbell (Reference 17), re- 
vealed that the ditti nal -variation factor (U) Is somewhat 
variable. A value of 32$ is considered valid for dates 
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prior to Fehrjiary 1963, and from August 1 963, onward, 1)650 

20% variation is considered valid. Between these dates, 

R is made to decrease linearly. 

Although in these equations the exponents m and 
n, which determine the mode of the longitudinal and 
latitudinal temperature variations respectively, arc 
kept distinct, it was found in practice that m = n. 

These values are not really known accurately and could 
be as small as 2.0. 

* 

. The constant B dptermincs the lag' of the tempera- 
ture maxi nun with respect to the uppermost point of the sun; . 
p introduces an asymmetry in the temperature curve whose 
location is determined by Y* 


Geomagnetic Activity 


To the temperature, T, which is calculated above, 
a correction must be added which accounts for atmos- 
pheric heating related to changes in the Earth's mag- 
netic field. The heating probably occurs in the E layer 
of -the ionosphere, but the mechanism involved is not 
well understood. The temperature correction, AT, is 
given by Jacchia, Slowey,. and Campbell (Reference 17): 


A T » 1.0° a p + 100* [l-cxpf-0.0$a r )3 


' u' 1 '!' 'u - 

(A 


where 





is ' the three -hourl y .planetary * geopaftrictic index/. 
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The quantity is a 
the earth's magnetic field 


measure of the variation in 
in a given three hour period. 


During magnetic storms the t«mperaturc changes 
generally lag behind the variations in a p by about five 
hours, due to conduction. There is some evidence of 
la.igcr temperature changes for given values of as 
one proceeds to higher geomagnetic latitudes. However, 
the amount of data indicating this is negligible at 
this time. 


The D650 subroutine allows for the magnetic 
heating effects with one modification. To minimize the 
input data for GEODYN, the 3 -hourly index (ap) is 
replaced' by a 24-! nurly or daily index (A^) . 

Generally, magnetic storms last for 2 or 3- days so 


that the temper 
a daily change, 

occur with a . 

.P 


a lure- calculation using !v trill reflect 

P 

but not the 3-hourly fluctations which 


The quantity A and the solar flux data is 
available from E.S.S.A., Boulder, Colorado. The pub 3 i 
cation is, "Solar Geophysical Data, Part I," 


D650 


Accurate daily values for both the solar and 
geomagnetic flux arc required for the computation of 
the exospheric temperature. In GEODYX, these values 
arc input via a BLOCK DATA routine, I KPT. This infor- 
mation may ;>c updated (cf subroutine ADEU1X) using the 
appropriate GEODYN Input Cards, The user should be aware 
of the fact that these tables are expanded os new inform, ' 
nation becomes available.. 


FlUXS 
FLUX A? 
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At the beginning of each run, a file is generated JANTIN'. 
for each satellite arc which contains the required flux 
data for the tine span indicated. Subroutine. JANTI1G 
is the routine v;hic.h sets up the flux tables, including 
averaging the daily values of solar flux over two solar 
rotation periods. The reason for this is to free the 
large amount of computer storage required for daily 
flux values over six years. As a matter of reference, 

.the associated COMMON BLOCK is PRIORI. 

8. 7. 2. 3 Tjio Density Computation 

• *. 

The density computation in GEODYN subroutine D650 

D650 is based on the density distribution versus 

altitude and exospheric temperature presented in Table 12, 

'■which is reproduced from Jacchia's . 1905 paper (Reference 
This data was obtained by numerical integration of 

the diffusion equation using an empirical temperature 
profile for each indicated exospheric temperature. 

This vast quantity of information was fitted 
(by v;OLF) to various degree polynomials of the form: 


LOG 


10 P D 


■E 

i 


y a. • T 
Z~y ij 


CM) h (i-l) 


( 7 ) 


where 


p0 *is the density, 






r 


04-48 
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T is the exospheric temperature, D6S0 

h is the spheroid height of the satellite 
(altitude), and 

a is a set of appropriate coefficients 



Unfortunately , a single polynomial of the type 
•presented is not completely descriptive. An examina- 
tion of Table 1 reveals that density is nearly in- 
dependent of temperature for low altitudes, but 
becomes increasingly dependent for height's above 
160 ton. Accordingly, appropriate' polynomials were 
chosen to account for the varying dependency of the 
variables. This necessitated the separation of 
Table 1 into three parts . 

The lower region (120 km - 360 kn) is expressed 

as a second degree polynomial which is solely a function 

of altitude. This is due to the fact that density is 

* 

not appreciably dependent on tcr.pcral.ure, in this region. 
The remaining regions of 160 km to 420 ton and 42C km to 
1000 km arc described by polynomials of fourth degree 
in both temperature and altitude. 

The coefficients for the selected polynomials are 
presented in Table 6. Theso coefficients have been 
modified to compute the natural log rather than tho 
decimal log ex the density . p^M^tjaBEOT Off TEE 




August 11, 1973 


The densities produced by these fitted polynomials 
differ from the densities in Table 1 by an RMS of 3.7 
percent. However, the fit does vary in different regions 
of the tab: c. In the region of verst fit, where the 
temperature it relatively low (700-1000* K) and the 
altitude varies fron 620-840 km, the RMS is somewhat greater 
being about 0.5 percent. The largest percent difference 
between densities is 13.2 percent and falls within the 
region described. 


:r;o 


The fits above could be improved by either going 
to' higher degree polynomials or by additional segmenta- 
tion of the table. However, these fits arc considered 

to be as accurate as the model being used. 

*■ 

v m * 

For satellite altitudes above 1000 km, the density 
is computed according to the extrapolation formula «iven 
by Jacchia (Reference 12) : 


Pn B P« + (Pjoo-pJc 


(b(h-1000) ] 


( 8 ) 


where 


Pi 8 e s 


REPROD UCIBlUTy OF THE 
ORIGINAL PAGE IS POOR 

air <*» Pp) as evaluated at 1000 km. 


is <i limiting value for the density . 
This is aero in subroutineORKSTY. 
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8. 7. 2. 4 Density Partial Derivatives DOS') 

VJSVAI. 

In addition to the density, GEODYN also requires 
the partial derivatives of the density with respect 
to the Cartesian position coordinates. These partials 
arc used in computing the drag contribution to the 
variational equations. 

As demonstrated above, the density is given by 


2 x * 

P D e exp (C 0 + c x h + C 2 h + C 3 h j (1) 

* ». • 

where ~ • 

# 

h is the spheroid height, and the 
C’j are coefficients which are polynomial s in 
temperature . 


We then have 


REPRODUCIBILITY OP THE 
original: PAGE IS POOR 


& Pn ./ 3h 

— * P n (C, + 2 C. h + 5 C t V) — 
8r u 1 * * 


(-J 


where 


t\ is the; t.m .of 
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The parlial derivatives —— arc computed in subroutine 

VEVAL. The quantities h, Pj } , and the arc computed 
in D5SO and passed through COMMON CLOCK DRCBLh. 
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8.8 


TIDAL POTENTIAL 


TIDAL 


The gravitational potential originating from solid 
earth tides caused by a single disturbing body is given 
(Reference 11). 



and the resultant acceleration on a satellite due to this 
potential is 


vu, 


*2 CM,. R e 


where 



(R d • h 2 ] r ♦ 6(R d • r)R d j 


( 2 ) 


kj is the tidal coefficient of degree 2 called the 
•"Love Number." 

G is the universal gravitational constant 


M_ is the mass of the earth. 


R is the moan earth radius. 

9 




M 4 




TO > 


Mj is the muss of the disturbing body. 


TI1HL 


M is the mass of the earth, 
e 

R, is the distance from COM* to 
d e 

r is the distance from CONI to 

e 

, i»s the unit vector frem COM^ 
G e 

A 

v is the unit vector from COM^ 

e 


COM** . 

vhe satellite, 
to COM d . 
to satellite. 


* Center of mass of the earth. 

**Center of mass of the disturbing body 


SECTION 9.0 

INTEGRATION AND INTERPOLATION 


GEODYN uses Cowell's Sum method for direct 
numerical integration of both the equations of motion 
and the variational equations to obtain the position 
and velocity and the attendant variational partials 
at each observation time. The integrator output is 
not required at actual observation times; it is 
output on an even integration step. GEODYN uses an 
interpolation technique to obtain values at the 
actual observation time. The specific numerical, 
methods used in GEODYN for this integration and 
interpolation are presented below. These procedures 
are controlled by subroutine ORBIT. 


9.1 INTEGRATION 

Let us first consider the integration of the 

equations of motion. These equations are three 

second order differential equations in position, and 

may be formulated as six first order equations in 

position and velocity if a first order integration 

scheme were used for their solution. For reasons of 

increased accuracy and stability, the position vector 

r is obtained by a second order integration of the 
•• • 

accelerations r, whereas the velocity vector r is 
obtained as the solution of a first order system. 
These are both multi-step methods requiring at least 
one derivative evaluation on each step. 



COKE IX 


The integration scheme is equivalent to the inter- 
polator with arguments 1 and 0 for predictor and corrector 
respectively. 

To integrate the position components, the predictor 
q 

r n+l = (S 2 + £ Y p r n-p )h2 (1 

p=0 


is applied, followed by a Cowell corrector: 



q 


< S 2 *£ 

p*0 


* 



r n-p+l 



C2) 


The velocity components are integrated using the 
predictor; 


q+1 

ll 3 < s l * E »p V p )h 

p=0 

followed by an Adams -Moulton corrector; 

R + 1 

r , ■ (S t + S' 8 F . )h • (4) 

n+1 1 1 P p n-p+1' v J 

p a 0 

In these integration formulae, h is the integration step 

* * 

size, q has the value ORDER- 2, and Yp> v p* Bp and 8p are co- 
efficients whose values, are obtained from subroutine COWCOF.* 


* Published numbers are in Reference 1. 
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Under certain conditions, a reduced form of this 
solution is used. It can be seen from the variational 
and observation equations that if drag is not a factor 
and there are no range rate, doppler, or altimeter rate 
measurements, the .velocity variational parti a Is arc 
not used. There is then no need to integrate the 
velocity variational equations. This represents a 
significant time saving. In the integration algorithm, 
the B matrix is zero and (1-H) is reduced to a three 
by three.' * 

PRECEDING PAGE BLANK NOT FILMED 
A detailed description of the H matrix and the 
X n and V n vectors can bo found in pages 16, 17 of 
Reference 2. 

Backwards integration involves only a few simple 
modifications to these normal or forward integration 
procedures. These modifications are to negate the 
step size, and invert the time completion test. 

The above integration procedures are implemented 
in GEODYN in the subroutine COWELL. The inversions 
for backwards integration are performed by COWELL and 
ORBIT. The matrix inversion is performed by subroutine 
DNVERT. 

THe default step size for these integration pro- 
cedures is selected on the basis of perigee height and 
the eccentricity of the orbit. The default step size 
selection is explained in detail in the Operations Manual, 
Volume III of the GECDYN System Documentation. This may 
be reset to some other fixed value on input. (See the 
STEP control card description in the above manual.) 
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COWELL 

DNVERT 




Variable Step Mode 


There is an optional variable step mode which is 
the default mode for high eccentricity orbits. The 
selection of ll» is mode of operation, it' default initial 
step size, halving error bound, and doubling error bound, 
or variable increase or decrease of step size are also 
explained in Volume III with the STEP control card. 


In the .variable step mode, the local error is COWELL 

compared against upper and lower error bounds to determine REARG 
whether the step size should be increased or decreased. 

This local error is computed as the difference between the 
predicted and corrected values of position. Both'the 

increasing and decreasing procedures require the table: 
back values of accelerations to be modified so as to be REARG 

compatible with the new step size. The decreasing requires 
the interpolation for mid-points. This interpolation is 
of course. on the back position, velocity and acceleration 
values. The increasing is achieved by discarding every 
other time point in the table of back values and then the 
refinement using the decreasing algorithm. 

% 

It should be noted that 2(0RDER-1)-1 of back values 
are saved when GEODYN is operating in v Triable step mode. 
Increasing of the step size is disabled for the following 
, ORDER-2 steps after a step size change; i.e., until the 
table of back values is again filled. 

These increaving and decreasing procedures are con- 
tained in subroutine REARG. 
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9. 2 


Till: INTEGRATOR STARTING SCHEME 


The predictor-corrector combination employed to 
proceed with the main integration is not self -starting. 

That ji?j each step of the integration recmires the 

•A *-* * 

knowledge of past values of the solution that are not 
available at the start of the integration. The method 
presented here is that implemented in the GEODYN 

subroutine START. 

# ¥ 

A Taylor series approximation is used to predict 
initial values of position and velocity. With these start- 
ing values, the Sum array is evaluated using epoch positions 
and velocities. Now the loop is closed by interpolations 
for the positions and velocities not at epoch and tl . „r 
accelerations evaluated. The Sunn are new again evaluated, 
this proceduie continues intil the Sums converge to the 
desired accuraiy. 


9.3 INTERPOLATION 

GEODYN uses interpolation for two functions. The 
first is the interpolation of the orbit elements and var- 
iational par-ials to the observation times; the second 
is the interpolation for nud -points when the integrator is 

4 

decreasing the step size. 

The formulas use'* by INTRP are: 


START 


INTR? 

COEF 


for positions and 






SECTION 10.0 

THE STATISTICAL ESTIMATION SCHEME 


The basic problem in orbit determination is to 
calculate, from a given set cf observations of the 
spacecraft, a set of parameters specifying the 
trajectory of a spacecraft. Because there are gen- 
erally more obsei vations than parameters, the parame- 
ters are overdetermined. Therefore, a statistical 
estimation scheme is necessary to estimate the 
'’best" set of parameters. 

The estimation scheme selected for GEODYN 
is a partitioned Bayesian least squares method. 

The comp, ete development of this procedure is pre- 
sented in this section. 

It should be noted that the functional re- 
lationships between the observations and parameters 
are in general non-linear; thus an iterative, pro- 
cedure is necessary to solve the resultant non-linear 
normal equations. The Ncwtcr-Raphson iteration 
formula is used to solve thcr.o equations. 
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10.1 


BAYESIAN LEAST SQUARES ESTIMATION'* 


Consider a vector of N independent observations 
z whose values can be expressed as known functions of 
M parameters denoted by the vector x. The following 
non-linear regression equation holds: 


z * f (x) ♦ o, (1) 


where o is the N vector denoting the noise on the, ob- 
servations. Given z, the functional form of f, and 
the statistical properties of o, we roust obtain the 
estimate of x that is "best” in some sense.** 

Bayes theorem in probability holds for proba- 
bility density functions and can be written as follows: 


P(x|z) 


P(x) 

p(z) 


p(z|x) . 


( 2 ) 


.«v & e 


P(«Sli.) is the joint conditional probability 
density function for the parameter vector x, given 
that the data vector z has occurred - 


^Vector notation in this sectic.i is that used by 
statisticians; i.e., an underscore denotes a vector. 

The symbol denotes the "best" estimate of the 
superscripted quantity. 

**For a complete discussion of the properties of estima- 
tors see Maurice G. Kendal? and Alan Stuart, Reference 1. 
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p(x) is the joint probability density function 

for the vector x; 

p(z) is the joint probability density function 

for the vector z; 

and 

P :z|x) is the joint conditional density function 
for the vector given that x has occurred; 

p(x) is often referred to as the a priori density 
function of x, and p(x|z) is referred to as the 
a posteriori conditional density function. In 
any Bayesian estimation scheme, we must determine 
this a posteriori density function and from this 
function determine a •'best" estimate of x, which 

a ““ 

can be denoted x. 

To obtain the a posteriori conditional density 
function, we must make an assumption concerning the 
statistical properties of the noise on the observations: 
the noise vector o has a joint normal distribution with 
mean vector £ and a variance -covariance matrix 
^ z is an NxN matrix and is assumed diagonal, that is, 
the observations are considered to be independent and 
uncorTelated. The "best" estimate of x, x, is defined 
as that vector maximizing the a posteriori density 
function; this is equivalent to choosing the mean value 
of this distribution. An estimator of this type has 
been referred to as the maximum likelihood estimate in 
the Bayesian sense. (Reference 2) 
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A further assumption is that the a priori density 
function p(x) is a joint normal distribution and is 
wiitten as follows: 



where 


is the a priori estimate of he parameter 
vector. 


/ ,, is the a priori variance-covariance matrix 
associated with the a priori parameter vector. 
E a is an MxM matrix, which may or may not be 
diagonal . 

The conditional density function *jx) can be 
written as follows: 


P(2|£> r 


Dctfr; 1 ) 

— 7^ 

2u 


1 1 

2 


exp 


1 

2 


i-t (x) 


IT 


rim 


( 4 ) 


It can be shown that maximizing the a poster iori density 
function p(x{j:) is equivalent to maximizing the product 
P(x)p(il£) because the density function p(t) is a con- 
stant valued function. Further, this reduces to mini- 
mizing the following quadratic form: 

* a 
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This results in the following set of M non-linear 
equations: 

® T £ * Z (iA-s) * 0 (6 > 

where B is an NxM matrix with elements 


B 


NM 


3£n(^> 

8x.. 


. n «r\TT TtTTtT 


This equation defines the Bayesian ieait squares 
estimation procedure. Ke have not . slated how the 
a priori parameter vector and vari anc/j- covariance 
matrix were obtained. In practice these a priori 
values arc almost always estimates that have been 
obtained from some previous data. In these cases 
the Bayesian estimates are identical to the classical 
maximum likelihood estimates that would be obtained 
if all the data were used; in this context the 
a priori parameters can be considered as additional 
observations . 




w 


i 



/ 



! 


C 


I 

) 

) 

I 
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The variance -covari ance matrix of x, V, is 
given by the following formula: 


V 



-1 


( 7 ) 


Solution of the Estimation Formula 


Equation 6 defines a set' of M non-linear equa- 
tions in M unknowns x*» these equations are solve'd 
using the Nevton-Raphson iteration formula. Equation 6 
can be written as follows: 




The iteration formula is 


f* 




i(nn) . 





( 8 ) 




rvhere 

x ^ is the n ** 1 approximation to the true 
solution x. 
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F(?) = B' 


E l ~ \ C“"' / \ . 

i-£(s) * L. x A -s) 

z v f A \ ' 


0 (9) 


Then diffe v entiating and neglecting second deriva- 
tives , 

REPRODUCES ILIT i OF VM 

r "I ORIGINAL PAGE IS POOR 

(¥)• K4 • £ : » 


(10) 


Substituting equation 10 in equation 8 gives 


i("U) . *(n) . (b^'b *rT fiT E 1 (i-it2) Cn) ) 


E t - ;« 


(ii) 


Now let x( n+1 )-x (n \ the correction to the n^Happroxi- 
mation, bo denoted by dx^ n+1 \ and let z-f(x^), the 
vector of residuals from the n " approximation, be 
dz ^ . Equation 11 becomes 


-1\-1 1 

*1 

-i 

/ 

£ 1 

A / | 

.b t 2) 

\r> / 
L~j ( 


Z 

* * 
/V 

(12) 
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10.2 


Till: PARTITION Li) SOLUTION 


In a multi-satellite, multi-arc estimation program 
such as GLOJYN, it is necessary to formulate the estima- 
tion scheme in a manner such that the information for 
all satellite arcs are net in core simultaneously. The 
procedure used in GEODYX is a partitioned Bayesian 
Estimation Scheme which requires only common parameter 
information and the information for a single arc to be 
in core at any given time. The aevelopmcnt of the 
GEODYN solution is given here. 

The Bayesian estimation formula has been devel- 
oped in the previous section as 


where 


dx 


(nM) 


[h T KB + V A * T j 1 B T Kdm + V’ 1 


-A * £ 


C r 0] 


SSSBRW 

is the a priori estimate of x. 

is the a priori covariance matrix associated 
with x^. 


CD 


W is the weighting matrix associated with the 
observations . 


X< n) is the approximation to x. 


dm is the vector of residuals (0-C) from the 
approximation. 


ESTIM 
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dx 


, (n + 1) 


is the vector of corrections to the 
parameters; i.c.. 


ESTIM 


x” +1 



dx 


(n+1) 


6 


is the matrix of partial derivatives of the 
observations with respect to the parameters 
where the i, element is given by 



The iteration formula given by this equation . 
solves the non-linear normal equations formed by mini- 
mizing the sum of squares of the weighted residuals. 


We desire a solution wherein x is partitioned 
according to a, the vector of parameters associated 
only with individual arcs; and k, the vector of parame- 
ter^ common to all arcs. For geodetic parameter esti- 
mation a consists of the sets of orbital elements, 
satellite parameters, and measurement biases associated 
with each arc, whereas k consists of the geopotential 
coefficients and station coordinates. 

As a. result of this partitioning , we may write 
B, the matrix of partial derivatives of the observations, 
as 
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a 


where 


;:sti 


3m. 

T 



We may also write V^, the covariance matrix of 
the parameters as 



( 3 ) 


where we have assumed the independence of the a priori 
information on the arc parameters and common parameters 
(in practice valid to an extremely high degree) . 
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Wc may now rewrite our iteration formula: 


ESTIM 



WR + V 
a a 


KB,. 

IV 

,r 

i<B k 

i 

“k * 


B^ Wdm (n ' + V^(a - 

T (»0 \ 
B. Wdm + V. (K ln) - kj 


A A K T ] 

aI * c k 


The required matrix inversion is obtained by 
partitioning. We write 


N 1 »J • A A k 


I T 

N 2 N 4 


A T K 


and, solving the resulting equations, determine 


N 1 • A' 1 ♦ [A 1 A k ] N 4 [aJ A' 1 ] 
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(7) F.STIM 


N. 


-A 


A k N 4 


and 


N , 


K 


* A,. A 



( 8 ) 


There is no problem associated with inverting A 
because the existence of the a priori information- alone 
guarantees this property. On the other hand, the 

T -i 

inverse of K - Aj. A A^ is not guaranteed to exist. 
High correlations between the parameters could make 
the matrix near singular. In practice, however, the 
use of a reasonable amount of a priori information 
eliminates any inversion difficulties. 


The iteration formula may now be written as 


da 

dk 




( 9 ) 


or 


da 


A' 1 ♦ (A -1 V) n 4 ( A k A ’ X ) 


dk . -x 4 Aj a' 1 C n ♦ N’,| C k 


C a ’ A '* A k »« C k 
( 10 ) 
( 11 ) 
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Noting the similarities between dji nnu dk, we write 


ESTIM 


da = A' 1 C a - A' 1 A k dk 


( 12 ) 


and rewrite dk as 


dk • N 4 (C k - AT A' 1 C a ) . 


( 13 ) 


Note that most of the elements of A aTe zero 
because the measurements in any given arc are inde- 
pendent - of the arc parameters of any other, arc. Also, 
the covariances between the a priori information 
associated with each arc is assumed to be zero. Thus 
both A and V & are composed of zeroes except for matrices, 
A r and V r , respectively, along the diagonal, where 


r is a subscript denoting the r** 1 arc, 
* 

• e.g., a r 








( 14 ) 


where £ ranges over the measurements in the r tft 
arc and i, j range over the parameters in the 
r th arc, a f . 
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ESTI. 


V is the partition of V_ associated with the 
r -n, a 

r* arc. 

The reader should note that A like A, is composed of 
zeroes except for matrices A r * along the diagonal. 

We shall also require the partitions of A,, and 

K 

C according *to each arc. These partitions are given 

ft 

by 









(15) 



(16) 


where the subscript r again denotes the r^ 1 arc and t 
ranges over the measurement partials and residuals in 
the r** 1 arc. 


Let us now investigate the matrix partitions 
in the solutions for da and dk. Wo consider A~* to be 
a diagonal matrix with diagonal elements A ** and C a 
to bo a colur. n vector with elements C r * Hence 


ESTI‘1 


A' 1 C. 


= rt C 
r r 


( 17 ) 


Jr 


■ b 


is the r , ' , ‘ element of the product matrix. Aj. is con- 
sidered to be a column vector .. ith elements A f j., thus 


» T .-1 
A k A '"a 


A rk A r^ C r 


(18) 


_ •% 

The elements in the product A x Aj. are given by 


-1 


L 


* V A rk 


kuwduciwutv op THE (19) 
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T -1 

We also require the product A A^. Its elements 
are given by 


T* - 1 

4 A A k 


A T a -1 a 
A rk A r A rk 


( 20 ) 


•J r,r 


The solutions for da and dk may now be rewritten 
taking into account the partitioning by arc: 



(21) 


dk 



v > 

jLj 


T ,-l 
rk ‘ 



( 22 ) 
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h'hc . c 


ESTIM 




( 23 ) 


These solutions form the partitioned Bayesian estima- 
tion scheme used in GLOOYN’. 

Additionally, the covariance matrix for the arc 
parameters must be updated to account for the simultan- 
eous adjustment of the common parameters: 


W ■ a;1 4 h 1 ** ( 


T 

a . 


( 24 ) 


Summary 


REPRODUCIBILITY OB’ THE 
ORIGINAL PAGE IS POOR 


The procedure for computer implementation is 

illustrated in Figure 1. This procedure is: 

# 

1* Integrate through each arc forming the 
matrices A r , A^, and C r ; and simultan- 
eously accumulate into the common 
parameter matrices K and Cj,. 

2. At the end of each arc, form 


C r 


( 25 ) 


and modify t.he common parameter matrices 
as fellows: 
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K ' ■ * ■ A rk \ l A rk 


HSTIM 


( 26 ) 


C k " 


( 27 ) 


The matrices <2a£, A r jj» and \ must also 
be put in external storage. 

After processing all of the arcs; i.e., • 
at the end oi a global or "outer" iteration, 
determine dk. Note that K has become 
and has been modified so that 


K" 1 C, 


( 23 ) 


The updated values for the common parameters 
are of course "iven by 


l-CR’-l) „ i-C”) ♦ at 


( 29 ) 


The arc parameters arc then updated to 
account' fw- the simultaneous solution of 
the coir.: o\ x ' ^meters. Information for each 
arc i? th\'.u\ in term; that is, the previously 
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- ■*> #orr«r - 


stored da 1 A , , and A The correction 
— r rk 

vector to the updated arc parameters is 
given by 


F.STIM 


da r * da; - (A; 1 A rk ) dk 


and hence 


a r (n+1) - ir Cn) ♦ da r 


( 30 ) 


(31) 


The covariance matrix for the aTC parameters, 
A; 1 , is updated oy 

a; 1 = a ; 1 * (A r x A rk ) K - 1 (A rk a; 1 ) (32) 


This completes the global iteration. 

It should be noted that if only the arc parameters 
are being determined, as is the case for "inner’' itera- 
tions, the solution vector is da| and hence the updated 
arc parameters are computed by 


a (n+ 1 ) 
-r 


a + da 1 
~r — r 


(S3) 
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Figure 1: Partitioned Intimation Procedure (Cont.) 
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Figure 1: Partitioned listimation Procedure (Cent.) 
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Figure 1: Partitioned Estimation Procedure (Cent.) 
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Fi^ rc 1: Partitioned Ur ti mat- ion Procedure (Cor.t.) 
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The common parameter matrix K is carried as a 
symmetric matrix. It is core-resident throughout the 
estimation procedure. Its dimension is set by the 
number of common parameters being determined and remains 
constant throughout tnc procedure. 


The arc parameter matrices A are also carried as 
symmetric matrices. Their dimensions vary from arc to 
arc according to the number of arc parameters being deter- 
mined. Only one arc parameter matrix A f and the^ corres- 
ponding covariance matrix A r j, are resident in core at 

any given time. These arc parameter matrices are stored 
on disk during step 2 of the above summary and recovered 
during step 3. 


The a p7ior i covariance matrix Vj. is not carried 
as a full matrix. The correlation c efficients be- 
tween each coordinate of a given station position are 
carried. The position coordinates of different stations 
and the g'copotential coefficients are considered to be 
uncorrclated . 


The a priori covariance matrices V are also not 
carried as full matrices. The drag coefficient, radi- 
ation pressure coefficient, and each bias arc considered 
to be uncorrelated. The covariance matrix for the epoch 
elements is carried. 


Tr-rulLlTY 
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liSTIM 


In terms of a subrout ine breakdown within GEODYN , 
this entire section is implemented in subroutine F.STIM 
with the exception of the matrix inversions. These 
inversions are done by s routine SYMIXV. 

0 * 

10.3 DATA EDITING 

The data editing procedures for GEODYN have two 

forms : 

• hand editing using input cards to delete 
specific points or sets of points, and 

• automatic editing depending on the weighted 
residual as component to a given rejection 
level. 

The hand editing is a simple matching of the GEOSRD 

appropriate GEODYN control card information with the DODSRD 

set of observations. This calling procedure is done 
in GLUDfN subroutines GEOSRD or DODSRD. 

The automatic editing of bad observations from NONAME 

a set of data during a data reduction run is performed 
in the GEODYN main program. Observations are rejected 
when 



t 


1 
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where 


i ^ ir 

WU.\.‘v.lC 

0 is the observation 

C is the computed observation 

o is the a priori standard deviation 

associated with the observation (input) 

k is the rejection level. 


The rejection level can apply either for all 
observations of a given type or for all observations 
of a given type from a particular station. This re- 
jection level is computed from 


k 


c 


M 


'R 


( 2 ) 


where 

is an input multiplier, and 

is the weighted RMS of the previous "outer" 
or global iteration. The initial value of 
Ev. is s 2 1 on input. 


It should be noted that both E M and E„ have default values, 



10.4 Electronic Bias 


t r* 


bsco:-.i 


For certain types of electronic tracking data 
(e.g., Dopptcr data), biases exist which are diffeier.t 
from one pass to the next. In many cases, these biases 
arc cf no interest per sc, although their existence 
must be appropriately accounted for 'if the data is to bo- 
used in an efrbit or geodetic parameter estimation. In 
addition, a single data reduction may have hundreds of 
passes of such electronic data, and the complete solution 
for each bias would require the use of an excessively 
large amount of computer core for storing the normal 
matrix for the complete set of adjusted parameters. 


The effects of electronic biases can be removed, 
with the use of only a small amount of additional core, 
based on the partitioning of the biases from the other 
parameters being adjusted in the Bayesian least squares 
estimation. The form which this partitioning takes can 
be seen from the solution of the basic measurement 
equation 


where 

6m * 

B Ab 
e 

■* BAX ♦ e Cl) 

sssssw* 


6m 

*s 

the vector of residuals (0-C) , 


Ab 

s 

ti e set of corrections that should be 
r.ac-7 to the cl. -tronic biases, 


B o 

31 

the matrix of partial derivatives of the 
mcr.sn’. encnls with respect to the biases. 
The elements of this matrix arc cither 


• 


l's or 0's. 
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RS 


i'.:; - the sot of corrections to be made to all 

the other adjustable parameters, 

B *” tne ci L ■ a u r pax lx al uciivativcs of the 

measurements with respect to die x 
parameters , 

£ c t.ic measurement noise vector. 


The least squares solucion o f (1) is 


A 

ftb 
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r T 
B e 

T 

B A \U 
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B T W6m "1 
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A 

Ax 


B T WB 

p 

T 

B W3 


B T N 6m | 
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1 w 

L 
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1 T 

with W the weight matrix (W = E(ee )), taken to be completely 

A 

diagonal in CEGDYN. The AX part of (2) can be shown to be 


h - [P T K® - Blv;E e 'JW 1 


(3) 


Jvo 


T 


T 


x [B W<Sm - 3 KB (B *KB A ) * B *W6m] 

C w 6 C 


To effectively remove the electronic bias effects, Eqn ' (3) 

states that the normal matrix must have B^WR (B ) * 

t c c 

T T 

B WB subtracted from it and the vector R W 6m must have 
c 

B T WB (B T WB )~ * Bottom subtracted from it. Due to the assumed 
e e e e 

independence of different measurements, it follows that these 
quantities which must be subtracted arc a sum of contribitfons 
for different passes, 
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where n^ is ‘the total number of passes with electronic 
biases and the subscript p denotes an array for measure- 
ments of pass p. The computation of the right hand sides 
of (4) and (S) requires the arrays 


T 

B„ N B « na x 1 array 
P P P 

B e T Ve f “ 1x1 a * r? ‘ y (&) 

T 

B^ h' = lxl array 

e P p p 

% 

where na is the numher of adjusted parameters other than 

biases affecting the arc in which the baises occur. Thus, 

• « 

na + 2 storage locations must be assigned for every bias 
which exists at any one time. 

The individual biases may be adjusted, based on the 
previous iteration orbital elements and force model para- 
meters. This bias can then be used, along with the above 
accumulated arrays to properly correct the sum of weighted 
squared residuals upon which the program does dynamic editing. 
Otherwise, however, it will not be possible for tlve statis- 
tical summaries to incorporate the adjusted values of the 

electronic biases unless substantial additional core is allocated. 
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SECTION 11.0 

GENERAL IXPUT/OilTJ’UT DISCUSSION 

GEODYN is a powerful yet flexible tool for 
investigating the problems of satellite geodesy and 

orbit analysis. This sane power and flexibility 

0 % 

causes extreme variation in both input and ouvput 
requirements. Consequently, GEODYN’ contains a great 
deal of programming associated with input and output 


11.1 INPUT 

There are two major functions associated with 
the input structure: 

These are the input of 

• Observation data, and 

% 

• GEODYN Input Cards. 

The observation data utilized by GEODYN in- 
cludes data from all the major satellite tracking 
networks. The observational types used to date, 
together with their originating networks and instru- 
ment types, sre: 

• Hight Ascension and Declination 

SAO Baker-Nunn cameras 

STADAN MOTS -cameras 



USAF 

PC- 1000 cameras 

USC5GS 

PX-4 cameras 

SPEOPT 

All of above except Baker -Nunn 


cameras 


o Range 


STADAN 

GRARR S-Band 
GSFC Laser 


SAO 

Laser 


AMS 

SECOR 

* 

C-Band 

FPQ-6 Radar 
FPS-16 Radar 


MSFN 

S-Band Radar 



• . Range Rate 

STADAN GRARR S-Band 

MSFN S-Band Radar 

•' Frequency Shift 

TRANET - Doppler 

a Direction Cosines 

STADAN Minitrack interferometer 

• X and Y Angles 

STADAN GRARR 

MSFN S*Band Radars 
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S' 


Azimuth and Elevation Angles 


STADAN 

C-BAND 


GS! ; C Laser 
FPQ-6 Radar 
FPS-16 Radar 


Time Delay and Fringe Rate 
C-BAND VLB I Radars 


; ) 
! I 


J 


; » 


c 


The Observations are required to be in either 
the format specified by the National Space Science 
Data Center (NSSDC) or the GSFC DODS System. 

The NSSDC format includes indicators to identify 
observation type, instrumentation source, reduction 
method, coordinate system, and information concerning 
tropospheric and ionospheric refraction corrections. 
Data in this format is input via subroutine GEOSRD. 

The DODS format includes indicators to identify 
observation type, satellite identification, ambiguity 
corrections, transponder channel when applicable, 
timing correction, and time reference system informa- 
tion. It also contains flags to indicate the need 
for transit timo correction or other types of pre- 
processing corrections. Data in this format is input, 
via subroutines DODSRD and DATBSE. 


GEOSRD 


DODSRD 

DATBSE 


H 


C 


The GEODYN Control Cards are the complete 
specifications for the problem to be solved including 
special output requests. Their input, controlled 
through subroutines ADFLUX and INOUPT, consists of 
data and perhaps variances for 

e Cartesian orbital elements 

• Satellite drag coefficients 
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ADFLUX 

INOUPT 
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• Satellite emissivity 

• Zero set measurement biases to be adjusted 

• Station positions 

o Geopotential coefficients 

• Surface densities 

• Earth tidal parameters 
and data for 

• Satellite cross-sectional area 

• Satellite mass 

o Integration times for the orbit 

• Epoch time of elements 

• Criteria for iteration convergence and data 
editing 

• Solar and geomagnetic flux 

Subroutine ADFLUX modifies the program internal 
data tables of solar and magnetic flux according to the 
input, requests. It also generates the scratch file of 
flux information to be used with each arc. 

Subroutine IK’OUPT interprets the (, EODYN Control 
Cords and sets the appropriate run parameters. It 

also generates the CEODYX run description and the 

* 

descriptions for all arcs. . \ 
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ADFLUX 

TN'OUPT 1 
1 



IN0L1PT 


Subroutine IN’OUPT references other routines to 
set up certain run parameters or to list selected run 
parameters in a particular format. 


It should be noted that the starting orbital . DODELM 

elements for some arcs may be recovered from the DODS 

Data Base by subroutine DODELM. 

* * 

11 . 2 Output 


Most of the output from GEODYN, not counting the NONAME 
descriptions of the input or run parameters, is pro- 0RB1 

duced by the NONAME program. Exceptions to this are SUMMARY 

the 0RB1 tape output, the residual summary and the run TYPOBB 

summary page. 


The printed output consists of a measurement 
and residual printout, residual summaries, and solution 
summaries as detailed below. 

For each arc: 

<• 

Measurement and Residual Printout 

• Measurement date 

• Measurement station 

• Measurement type 

• Measurement value 
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Measurement residual 


• Ratio to sigma 

• Satellite elevation 

Residual Summary by Station and Type 
c Station 

• Measurement type 

• Number of measurements 

• Mean of residuals 

• Randomness measure 

• Residual RMS about zero 

« Number of weighted residuals 

•. Mean ratio to sigma for weighted residuals. 

9 Randomness measure for weighted residuals 

• RMS about zero for weighted residuals 
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Residual Summary by Type 


• Measurement type 

o Numbe* of weighted residuals 

• Weigh tec! RMS about zero 

© Weighted i'MS about zero for all types together 
Element Summary 

• a priori Cartesian elements 

• Previous Cartesian elements 

• Adjusted Cartesian elements 

• Adjustment to Cartesian elements (delta) 

• Standard deviations of fit (sigmas) 

•. Position RMS 

• Velocity RMS 

• a priori Kepler elements 

• Previous Kepler elements 
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• Adjusted Kepler elements 

• Adjustment to Kepler elements (delta) 

• Double precision adjusted Cartesian elements 
(current best elements for arc) 

# t 

Adjusted Force Model Parameter Summary for Arc 

• Drag Coefficients, Solar Radiation Pressure Coef- 
ficient, and/or resonant geopotential* coefficients 

• a priori coefficient value 

• Adjusted coefficient value 

• a priori standard deviations for coefficient 

• Standard deviation of fit for coefficient 
% 

Adjusted. Parameter Summary 

• Instrument biases - timing bias and/or 
constant bias 

• a priori bias value 

• Adjusted bias value 

• a priori standard deviation for bias 
o Standard deviation of fit for bias 
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Time period of coverage 


The following items are printed on the last inner iteration 
of every outer iteration. 

o Apogee and perigee heights 

* fc 

• Node rate and perigee rate 

• Period of the orbit 

• Drag rate and period decrement if drag is 

being applied 

• Updated covariance matrix for. Cartesian arc 
elements 

• Adjusted arc parameter correlation coefficients 

After all arcs : 

Total Residual Summary . 

* 

t Total number of weighted measurements for 
each measurement type 

• Total weighted RMS for each measurement type 

• Total number of weighted measurements 

• Total weighted RMS 
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Station Summary 


• Earth-fixed rectangular coordinates and 

geodetic (t|>,X,h) coordinates 

• a priori coordinate values 

• * *a priori standard deviations for coordinate 

values 

• Adjusted coordinate values 

• Standard deviation of fit for coordinate values 


• Correlations between determined coordinate 
values 


Gecpotential Summary 

• C and S _ coefficients for each n.m set 

• nm nm ’ 

determined 

•' a priori values 

• Adjusted values 

• Ratios of adjusted sigma to a priori sigma 

for each coefficient 

• Standard deviations of fit for coefficients 
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Surface Density Summary 


* .* 
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• Surface Density Block Centers 

• Block Areas 

o a priori values 

• ' ‘adjusted values 

• a priori uncertainties 

* 

• adjusted uncertainties 

Arc Summary for Outer Iteration - For each arc 
a Updated Cartesian elements for arc 

• Correlation coefficients between individual 

arc parameters 

% 

t Standard deviation of fit for arc parameters 

• Correlation coefficients between individual 

arc parameters and parameters common to all 
arcs 


Common Parameter Correlation Coefficients 

• Geopotential coefficients 

• Cartesian station positions 

• Surface Densities 


U. 2-7 






) 


f 


> 


; I 

I > 

I 

l 

i i 


i 


I t 


c 


Ss 



GF.ODYN also produces an XYZ and Ground Track 
listing upon request. This is the normal printout for 
1- orbit Generation Mode. In addition an osculating ele- 
ment printout is provided on option. 

The tape output from GEODYN consists of 
s the 0RB1 tape, 

• the XYZ and Ground Track tape, 

• a DODS formatted data tape, 

t a binary residual tape 

© a simulation data tape. 

The XYZ and Ground Track tape and the binary residual 
tapes are used as input to GEODYN support programs. 
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11.3 Computations for Residual Summary 
% 

The residual summary information is computed in 
subroutine STAINF for printing by the main program. 

The formulas used in this subroutine for computing 
each statistic are presented below. 


STAINF 


The mean is 



11.3-1 



where 


/ 
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R. are the residuals 

n is the number of residuals 

is the number of electronic biases 

N. are the residuals contributing to the bias 
J computation 

b_. is the value of the electronic bias. 

v • 


( 
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The RMS is the square root of the sample 


variance : 


STAIXF 


where 


n. 

, r b 

- 1 E^.-Ev* 

» Lw j-i 1 


The expected value of the sample variance differs from 

2 

the population variance a : 


E(s 2 ) 


o - var (P c ) 


or rather 


E(s 2 ) 


° 2 u - b- 


Hence we may make a better estimate of o^’by computing 
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H This is known as Bessel's correction. This computed 
value for the standard deviation, o, is also called 
the RMS about zero. 


The randomness measure used in GEODYN is from a 
mean square successive difference test. V;e have 



( 6 ) 


when 


RND is the random normal deviate, our statistic; 
2 

s is the unbiased sample variance; anc 



1 

2(n-l) 


n-1 

E <*i.i 


i*l 





V 


V.. 


c 


2 

Note that d is the mean square successive difference.' 


For each i the difference R 


i+1 


Rj has mean zero 


and variance 2c under the null hypothesis that 

(R, ,....R_) is a random sample from a population with 

* ’ n 2 2 2 
variance o . The expected value of d is then o . 

If a trend is present d* is not altered nearly so much 

as the variance estimate s^, which increases greatly. 

Thus the critical region RNO constant is employed in 

testing against the alternative of a trend, (Reference 1) 
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STAIN 



STAIXF 


In order to use this test, of course, it is 
necessary to know the distribution of the RND. It 
can be shewn that in the case of a normal population 
the expected value is given by 


E (RND) - 1, 

0 * 


( 7 ) 


the variance is given by 


var (RND) 



( 8 ) 


and that the test statistic, RND r is approximately 
normal for large samples (n > 20) . 
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11 . 4 Kepler Elements 

The Kepler elements describe the position of the 
satellite as referred to an ellipse inclined to the orbit 
plane. This is shown in Figures 1 and 2. The definitions 
of these elements arc: 

# fr 

a - semi-major axis of. the orbit 
e - eccentricity of the orbit 

1 - inclination of the orbit plane 
ft - longitude of the ascending node 

- argument of perigee 
M - mean anomaly 
E - eccentric anomaly 

2 - true anomaly 





Apogee height and perigee height are sometimes used 
in place of a and e to describe the shape of the orbit. As 
can be seen in Figure 1, the radius at perigee is a(i-e) 
and tnat at apogee is a(l+e). The heights are determined 
by subtracting the radius of the reference elipsoid at the 
given latitude from the spheroid height of the satellite, 

•The computations of these last are detailed in section 
5.1. 
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Conversion to Kepler Elements 


ELEM 


The computation of Kepler elements from the 

• • • 

Cartesian positions and velocities x,y,z,x,y,z is as 
follows : 

Compute the ‘angular momentum vector per unit mass: 


F = r x r 


( 1 ) 


where r is the position vector and r is the velocity 

2 — 

vector. Note that v » r • r. The inclination is 
given by 


i 



( 2 ) 


From the vis-viva or energy integral we have 


• • 

V* - CM , U) 

where G is the universal gravitational constant and-M 
is the mass of the primary about which the satellite is 



orbiting. Thus wo have 


SLBM 


2 _ v_ 
r ' GM 


2 -1 


Recalling Kepler's Third Law, 


h 2 = GM a (1-e 2 ), 


we determine 


f. .rf 

1 \aGM / * 


The longitude of the ascending node is also 
determined from the angular momentum vector: 


Q ■ tan 



The true anomaly, f, is computed next, 
integrating 


Note that in 


r x IT ■ GM r/y 


one arrives at 
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r - L'-l ‘ V - 


(9) ELEM 


F x IT = GM (r + e) 


where e is a constant of integration of magnitude equal 
to the eccentricity and pointing toward perihelion. 

Thus , 


r x e » re sin f 



( 10 ) 


or, performing a little algebra, 


sin f 


2 — — 
a (1-e ) r * r 

reh 


( 11 ) 


The cosine of the true anomaly comes from 


r . «iksX , 

I+ecosT 


( 12 ) 


that is 


cos £ * 


er 


1 

e 


( 13 ) 


The true anomaly is then 


£ » tan 


‘1 /sin j \ 

\coT? / 


(M) 
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At this point a decision must be made as to whether the 
orbit is a ellipse (l>e>0) or a hyperbola (l<c< to ). 

For an elliptic orbit, the eccentric anomaly is computed 
from the true anomaly: 


cos E = 


cos f + e 
1+e cos f ’ 

/T c x si, f 

1+e cos t 


E * tan 


-1 / sin E \ 
l cos E ) * 


I 1 c 


The mean anomaly is then computed from Kepler’s 
equation: 


M * E - c sin E. 

In the case of a hyperbolic orbit, we use an 
equation analogous to Kepler's equation by introducing 
F, in place of E. The eccentric anomaly is the same 
as above ; 


where 


f> A , -l, si nh 1* ^ 

F * tanh t-TSSiTT 5 


sinh F * 


cosh F 


fi^c 2 sin f 
1+c cos f 
cos f + e 
1 + e cos f 


ttnODUUUULnY OP 1HB 
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The mean anomaly is 


(13b) 

M = e sinli F - F 
whore F «* *n [sinli F + cosh F] 

F is computed by using the definition of sinh and cosh 

' * 'F 

sinh F = ^ — 

2 

cosh F = ~ — ^—5 — 

2 

(sinh F + cosh F) * (e^-e"^+e^+e~^) 


The central angle u is the angle between the satellite 
vector and a vector pointing toward the ascending node: 


cos u 


* x cos 0 * y sin ft 
.* r 


( 19 ) 


sin u 


(y cos n • x sin fl) cos i * z sin i 

r 


( 20 ) 



( 21 ) 


The argument of perigee is then 


PARTIAL DERIVATIVES Oi .EPLER ELEMENTS 
The partial dcrivrr •, -:f Kepler elements with 

• ft 4 

respect to . -,y , ?. ,x , / , • ioIIows: 

P.D.* of inclination; 

~ # i .« A[B(y.h z -a.h y )-y] 

\y i ■ ;A[B(z.h x -x.h z )+x] 
i c AB(x.h y -y.h x ) 
i - A[B(z.h -y.h )+y] 

dx y z 

i * A[B(x.h -z.h ) -x] 

3y z x 

ir- i * AB(y.h -x.h ) 

3z A 

where 

1 

A * 

h. sin- i 

cos i 
B - 

h 

P.D. of semi-major axis; 

3 V.2a 2 

— a - 

3V r A 
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partial derivatives 
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1ft 


where 

V * x,y,z. respectively. 

* 2 
3 V.2a* 

_ — , - -,- 

3V & GM 

where 

A » 

• • • • 

V * x,y,z, respectively. 

P.D. of eccentricity: 

e * C[x.D-| (y.h z -z.h y )] 

|y e - C[y.D-| (i.h x -i.h 2 )J 

^ e - C[z.D-| (x.hy-y h x )] 

K e = C [x. D’ (z.h -y.h z )] 
3x 7 

“r e ■ C[y.D’ -j (x.h z ."S h x )) 
3y 

K e - C[z.D’-| (y.h x -x.h )] 

dz 


I 

ji 

i 

i 

i 
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. 1 


where 

? 

I 


r - 1 

c " eTGM 
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P.D. of ascending node: 


CLEM 


C 


ox 

3_ 

3y 


ft = H. z . h 


ft = 11. 2. h 




^ ft 
ax 

a 


% • 


-ll(y.h y -x.h x ) 


* -11. z. h. 


3y 

a_^ 

az 


ft = -n z .h y 


ft ■ H(y.h y +x.h x ) 


vhere 


H 


1 

m 


h.. +h 

A 


P.D. of mean anomaly: 


3 M «* I S - sinE 3e 


37 


5v 


where 


i [ 


* 

i 


S * ( l ( 37 ’ I * cos E ||) / o sin E 


3e> 


v represents x, y, z. *, y, t respectively, end |i-|, 


o, 21 - o, II 
ay az 
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P.D 


of the argument of perigee 


EL 1 IM 


_i w = H - II 

3 v 3 v 3 v 


where 



9 $ (tan 


1 


sin u % 
cos u ' 


cos^u 


3u 

3 V 


sin u 


• i 

3 v 


cos u 


3 U 
3 V 


-1 

sm u 


3 __ 

3 V 


cos u 


and 


3 f a 

*v. 


3 


(tan 


1 sin f 
cos £ 


) 


simiiarily 


and 


li 


-1 


3V Fir. f ax 


co 3 f 


-g^cosu- ( (xs inft -ycosD) cosa|£ - s 


to*CQS 
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In GEODYX, this conversion from x,y,z,x,y,z to a,e,i, 
«o,M and the partial derivatives are performed by subroutine 
ELEM. 


Conversion From. Kepler Elements 

The input elements arc considered to be a,e,i, 
and M and the Cartesian elements are required. 

An iterative procedure , Newton’s method , is 
used to recover the eccentric anomaly. For an ellip- 
tic orbit, the iterative procedure is, from Kepler’s 
equation (M « E - e sin E) , 

E* * E - E ~ e sin E - M 
1 - e cos E 

For a hyperbolic orbit, the iterative procedure is 

p, m p e sinh F - F - M 
e cosh F - 1 

% 

where F, sinh F and cosh F are defined previously. 

This conversion procedure for converting 
a,e,i,n,w,M to x,y,z,x,y,i is performed in the GEODYN system 
by subroutine POSVEL. 


POSY EL 



The vectors X and IT are computed. A is a vector 
in the orbit plane directed toivard peri center with a 
magnitude equal to the semi-major axis of the orbit: 


X = a 


cos w cos ft - sin oj sin ft cos 
cos (o sin ft + sin w cos ft cos 
sin a) sin i 


i 

i 


5 is a vector in the orbit plane directed 90° counter 
clockwise from X with a magnitude equal to the semi- 
minor axis of the orbit. 


S. - a 


sin (i) 
sin u) 


cos ft - cos w sin ft cos 
sin ft + cos w cos ft cos 
cos a) sin i 


The position vector r is then 


r » (cos E - e) S + (sin E) S’ 


The velocity vector is given by 


r » E j(-sin E) 


X ♦ (cos E) F 


t 

] 


i POSVEL 
i (24) 


(25) 


( 26 ) 


where E is given by 
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Node Rate 


and Pen 


/t AO 


Rate 


The node rate Q and perigee rate w are computed from 
Lagrange’s Planetary Equation?. As these are for printout . 
only, GEODYN uses just the Earth oblateness term in the 
geopotential. From Reference 4, page 39, we have 


ft = 


(d 



COS 1 

(l-e 2 ) j 


(1-5 cos 2 i) 
(W 2 ? 


Cl) 


( 2 ) 


in radians per second, or rather 


ft. - -9.97 


r 


COS 1 


(3) 


-4.98 I — 



-3.5 (1 „ 5 co$ 2 i) 




(*) 


in degrees per day. Tho quantities used in the atvqvc equations 
are defined as:’ , • f.V 


is the semi -major axis of the Earth 



i 

I 


t 


J 

1 J 

i I 



GM is the product of the universal gravitational 
constant G and the mass of the Earth M 

C ?n is the Earth oblateness term in the geo- 
potential (see Section 8.3). 

a semi-major axis of the orbit 

'e eccentricity of the orbit 

i inclination of the orbit 


» 
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Period Decrement and Dr^g Rate 


< 

The period decrement and the drag rate are determined 
from the partial derivatives of the position and velocity 
with respect to the drag coefficient at the final integrator 
time step in the given arc. These (multiplied by the drag 
cocf f icicrit) ’ represent the sensitivity of the position or 
velocity to drag effects. Let us define 


( 1 ) 

where 

r is the satellite (inertial) position vector 
Cp is the drag coefficient 
We also define 



d 

E& * (r) • C 


The (two-body) period of the orbit is 



( 2 ) 


( 3 ) 


C 
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where 


Thus 


a is the semi -major axis of the orbit 

GM is the product of G, the universal gravitational 
constant, and M, the mass of the Earth. 


AP = 3 


■y 1 

V GM 


The vis viva or energy integral has 


GM r; ■ 


h^nce 


2 r • r 
r GM 


Recognizing that A(r) is AU and A(r) is SU, 


2 r - IC r • JU 

J p MW — M i ♦ , 

2 7 • ft r 3 GM 
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The effect of the drag on the period is then given by 



( 8 ) 


The daily rate or period decrement is computed as AP^ t 
where At is the elapsed time (in days) between the last 
integrator time point and epoch. 


The drag rate is computed from the along track 
'actually normal) portion of 7U, that is AD M , We need to 

A IN 

construct the unit vector along track, L. The velocity 
vector 7 may be resolved into a radial component and a 
component normal to the radius vector . The magnitude of 
the normal component is found by the Pythagorean Theorem: 


- ft • » (; r • r) 


(9) 


• A 

The unit normal vector L is then' 



The normal portion of AH Is then 


( 10 ) 



A 

AD^ * [• • 3Hj 


( 11 ) 



This AD^, represents the along- track position effect 
due to drag over the integrated time span. The drag rate 
is computed as AD^, A t 2 where At is again the elapsed time 
in days. 
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